// Author : Zhuowen Tu


// MCMC (Markov Chain Monte Carlo) region segmentation
#include "mcPublic.h"
#include "mcmcLayerRegionElements.h"
#include "mcmcLayerRegion.h"
#include "mcmcBlock.h"
#include "mcRaster.h"
#include "mcmcRegionImage.h"
#include "mcmcException.h"
#include "mcHeap.h"
#include "mcmcImageFilters.h"

#define	MINI_DIST1		3.0
#define	MINI_DIST2		2.0
#define	MINI_SIMP_DIST	0.6
#define MU				3.0
#define FORCE_NORMAL	5.0

void VerticesNearList::Insert(MCMC_Point &ver0, MCMC_Point &ver1)
{
	VerticesNear ver_near;
	ver_near.ver0 = ver0;
	ver_near.ver1 = ver1;

	list.Insert(ver_near);
}

bool VerticesNearList::Exist(MCMC_Point &ver0, MCMC_Point &ver1)
{
	bool			bValid;
	VerticesNear	ver_near;

	bValid = list.GetFirst(ver_near);
	while (bValid)
	{
		if ((ver_near.ver0 == ver0 && ver_near.ver1 == ver1) || (ver_near.ver0 == ver1 && ver_near.ver1 == ver0))
			return true;
		bValid = list.GetNext(ver_near);
	}
	return false;
}

///////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////
// manage nodes for vertices, contour segments, regions, and groups

/////////////////////////////////////////////////////////////
MCMC_MemoryNodeManagement::MCMC_MemoryNodeManagement()
{
	num_of_nodes	= 0;
	head_node_new	= 0;
	head			= -1;
	tail			= -1;
	current			= -1;
	max_space		= 0;
	p_nodes			= NULL;
}

MCMC_MemoryNodeManagement::~MCMC_MemoryNodeManagement()
{
	p_nodes	= NULL;
}

void MCMC_MemoryNodeManagement::Empty(void)
{
	num_of_nodes	= 0;
	head_node_new	= 0;
	head			= -1;
	tail			= -1;
	current			= -1;
	max_space		= 0;
	p_nodes			= NULL;
}

ofstream & operator>>(MCMC_MemoryNodeManagement &n, ofstream &strm)
{
	strm.write((char*)&(n.num_of_nodes), sizeof((n.num_of_nodes)));
	strm.write((char*)&(n.head_node_new), sizeof((n.head_node_new)));
	strm.write((char*)&(n.head), sizeof((n.head)));
	strm.write((char*)&(n.tail), sizeof((n.tail)));
	strm.write((char*)&(n.current), sizeof((n.current)));
	strm.write((char*)&(n.max_space), sizeof((n.max_space)));

	return strm;
}

ifstream & operator<<(MCMC_MemoryNodeManagement &n, ifstream &strm)
{
	strm.read((char*)&(n.num_of_nodes), sizeof((n.num_of_nodes)));
	strm.read((char*)&(n.head_node_new), sizeof((n.head_node_new)));
	strm.read((char*)&(n.head), sizeof((n.head)));
	strm.read((char*)&(n.tail), sizeof((n.tail)));
	strm.read((char*)&(n.current), sizeof((n.current)));
	strm.read((char*)&(n.max_space), sizeof((n.max_space)));

	return strm;
}

MCMC_MemoryNodeManagement& MCMC_MemoryNodeManagement::operator=(const MCMC_MemoryNodeManagement& a)
{
	if (this != &a)
	{
		num_of_nodes	= a.num_of_nodes;
		head_node_new	= a.head_node_new;
		head			= a.head;
		tail			= a.tail;
		current			= a.current;
		max_space		= a.max_space;
	}

	return *this;
}

// operations
void MCMC_MemoryNodeManagement::Set(MCMC_MemoryNode *p_nodes1, const int max_num, const char by_type1)
{
	max_space = max_num;
	by_type = by_type1;

	p_nodes = p_nodes1;
}


void MCMC_MemoryNodeManagement::Init(MCMC_MemoryNode *p_nodes1, const int max_num, const char by_type1)
{
	p_nodes		= p_nodes1;
	max_space	= max_num;
	by_type = by_type1;

	num_of_nodes	= 0;

	head_node_new	= 0;
	head			= -1;
	tail			= -1;
	current			= 0;

	int		i;

	if (p_nodes != NULL)
	{
		switch (by_type)
		{
		case 0: MCMC_Node *p_mnodes;
				p_mnodes = (MCMC_Node*)p_nodes;
				for (i=0; i<max_space; i++)
					p_mnodes[i].next = i+1;
				break;

		case 1:	MCMC_Vertex *p_vertices;

				p_vertices = (MCMC_Vertex*)p_nodes;
				for (i=0; i<max_space; i++)
					p_vertices[i].next = i+1;
				break;

		case 2:	MCMC_Contour_Segment *p_csegments;

				p_csegments = (MCMC_Contour_Segment*)p_nodes;
				for (i=0; i<max_space; i++)
					p_csegments[i].next = i+1;
				break;

		case 3: MCMC_AtomicRegion *p_atomic_regions;

				p_atomic_regions = (MCMC_AtomicRegion*)p_nodes;
				for (i=0; i<max_space; i++)
					p_atomic_regions[i].next = i+1;
				break;

		case 4: MCMC_GRegion *p_gregions;

				p_gregions = (MCMC_GRegion*)p_nodes;
				for (i=0; i<max_space; i++)
					p_gregions[i].next = i+1;
				break;
		}
	}
}


int MCMC_MemoryNodeManagement::New(void)
{
	int id_new;
	id_new = head_node_new;

	if (id_new >= max_space)
	{
		cout << "new" << endl;
		cout << "id_new: " << id_new << "      max_space: " << max_space << endl;
		XPT_OutOfBound  *pxpt;

		pxpt = new XPT_OutOfBound("Out of memory");
		throw pxpt;

		return -1;
	}


	MCMC_MemoryNode	*pmemory_node=GetNode(id_new);

	assert(pmemory_node != NULL);

	// delete it from the garbage can
	head_node_new = pmemory_node->next;
	pmemory_node->next = -1;
	num_of_nodes++;

	// put this region to the active region list
	if (head == -1)
	{
		head = id_new;
		tail = head;
		pmemory_node->last = -1;
	}
	else
	{
		MCMC_MemoryNode *pmemory_node_tail;

		pmemory_node_tail = GetNode(tail);
		assert(pmemory_node_tail != NULL);
		pmemory_node_tail->next = id_new;
		pmemory_node->last = tail;
		tail = id_new;
	}
	return id_new;
}


bool MCMC_MemoryNodeManagement::Delete(const int id)
{
	if (id<0 || id>=max_space)
		return false;

	MCMC_MemoryNode	*pmemory_node=GetNode(id);
	assert(pmemory_node != NULL);

	// is this vertex in the list?
	if (pmemory_node->last<0 && head!=id)
		return false;

	if (head == id)
	{
		head = pmemory_node->next;
		if (head < 0)
			tail = head;	// empty
	}
	else
	{
		MCMC_MemoryNode *pmemory_node1,*pmemory_node2;

		pmemory_node1 = GetNode(pmemory_node->last);
		pmemory_node2 = GetNode(pmemory_node->next);
		pmemory_node1->next = pmemory_node->next;
		if (pmemory_node2 != NULL)
			pmemory_node2->last = pmemory_node->last;
		else
			tail = pmemory_node->last;	// reach the tail
	}
	pmemory_node->last = -1;
	pmemory_node->next = head_node_new;
	head_node_new = id;
	num_of_nodes--;

	return true;
}


/////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// dynamic space

// construction & destruction
MCMC_DynamicSpace::MCMC_DynamicSpace()
{
	max_space = 0;
	p_data = NULL;
	num_active = 0;
}

MCMC_DynamicSpace::~MCMC_DynamicSpace()
{
	Free();
}

void MCMC_DynamicSpace::Free(void)
{
	if (p_data != NULL)
		delete []p_data;
	p_data = NULL;
	max_space = 0;
	num_active = 0;
}

void MCMC_DynamicSpace::Init(void)
{
	// free allocated space
	Free();

	// allocate memory
	max_space = 3;
	p_data = new int[max_space];
	assert(p_data != NULL);
}

// operators
MCMC_DynamicSpace& MCMC_DynamicSpace::operator=(const MCMC_DynamicSpace&a)
{
	if (this != &a)
	{
		Empty();

		if (a.max_space > 0)
		{
			GrowSpaceTo(a.max_space);
			for (int i=0; i<a.num_active; i++)
				p_data[i] = a.p_data[i];
		}

		num_active	= a.num_active;
		current		= a.current;
	}

	return *this;
}

ofstream & operator>>(MCMC_DynamicSpace &s, ofstream &strm)
{
	strm.write((char*)&(s.max_space), sizeof(s.max_space));
	strm.write((char*)&(s.num_active), sizeof(s.num_active));
	if (s.num_active > 0)
		strm.write((char*)(s.p_data), sizeof((s.p_data)[0])*s.num_active);

	return strm;
}

ifstream & operator<<(MCMC_DynamicSpace &s, ifstream &strm)
{
	int	temp;

	strm.read((char*)&temp, sizeof(temp));
	s.GrowSpaceTo(temp);

	strm.read((char*)&(s.num_active), sizeof(s.num_active));
	if (s.num_active > 0)
		strm.read((char*)(s.p_data), sizeof((s.p_data)[0])*s.num_active);

	return strm;
}

// operations
void MCMC_DynamicSpace::GrowSpaceTo(const int max_num)
{
	int *ptemp,i;

	if (max_space<max_num && max_num>0)
	{
		ptemp = new int[max_num];
		if (p_data != NULL)
		{
			for (i=0; i<num_active; i++)
				ptemp[i] = p_data[i];
			delete []p_data;
		}
		p_data = ptemp;
		max_space = max_num;
	}
}

bool MCMC_DynamicSpace::Exist(const int id) const
{
	int i;

	for (i=0; i<num_active; i++)
		if (p_data[i] == id)
			return true;
	return false;
}


bool MCMC_DynamicSpace::Add(const int id)
{
	if (id < 0 || Exist(id))
		return false;

	if (num_active >= max_space)
		GrowSpaceTo(max_space*2);

	assert(num_active < max_space);
	if (p_data == NULL)
		return false;

	p_data[num_active] = id;
	num_active++;

	return true;
}

bool MCMC_DynamicSpace::Delete(const int id, const int id_replace)
{
	int	i,j;
	bool	breplace_exists=false;

	if (id_replace > 0)
		breplace_exists = Exist(id_replace);

	for (j=0; j<num_active; j++)
	{
		if (p_data[j] == id)
		{
			if (id_replace >= 0 && !breplace_exists)
				p_data[j] = id_replace;
			else
			{
				for (i=j; i<num_active-1; i++)
					p_data[i] = p_data[i+1];
				num_active--;
			}
			return true;	// success
		}
	}
	return false;
}


int MCMC_DynamicSpace::GetFirst(void)
{
	current = 0;
	if (current < num_active)
		return p_data[current];
	else
		return -1;
}

int MCMC_DynamicSpace::GetNext(void)
{
	current++;
	if (current < num_active)
		return p_data[current];
	else
		return -1;
}


//////////////////////////////////////////////////////////////
// MCMC Node
/////////////////////////////////////////////////////////////
// operators
MCMC_Node & MCMC_Node::operator=(const MCMC_Node &a)
{
	if (this != &a)
	{
		x = a.x;
		y = a.y;
		link[0] = a.link[0];
		link[1] = a.link[1];
	}

	return *this;
}

// write to a stream
ofstream & operator>>(MCMC_Node &n, ofstream &strm)
{
	(MCMC_MemoryNode)n >> strm;

	strm.write((char*)&((n.link)[0]), sizeof((n.link)[0]));
	strm.write((char*)&((n.link)[1]), sizeof((n.link)[1]));
	strm.write((char*)&(n.x), sizeof(n.x));
	strm.write((char*)&(n.y), sizeof(n.y));

	return strm;
}

// read from a stream
ifstream & operator<<(MCMC_Node &n, ifstream &strm)
{
	(MCMC_MemoryNode)n << strm;

	strm.read((char*)&((n.link)[0]), sizeof((n.link)[0]));
	strm.read((char*)&((n.link)[1]), sizeof((n.link)[1]));
	strm.read((char*)&(n.x), sizeof(n.x));
	strm.read((char*)&(n.y), sizeof(n.y));

	return strm;
}

double MCMC_Node::TotalCodingBits(void)
{
	double dbytes;

	dbytes = 6.0/8.0;		// x,y

	return dbytes;
}

//////////////////////////////////////////////////////////////
// MCMC Vertex
/////////////////////////////////////////////////////////////
MCMC_Vertex::MCMC_Vertex()
{
	player_region = NULL;
	p_contour_segments = NULL;
}

MCMC_Vertex::~MCMC_Vertex()
{
	Empty();
	player_region = NULL;
}

// operator
MCMC_Vertex & MCMC_Vertex::operator=(const MCMC_Vertex &a)
{
	if (this != &a)
	{
		if (a.Valid())
		{
			GrowDegreeTo(a.Max_Num_Vertex_Degree);

			for (int i=0; i<a.Max_Num_Vertex_Degree; i++)
				p_contour_segments[i] = a.p_contour_segments[i];
		}
		degree			= a.degree;
		x				= a.x;
		y				= a.y;
		current_cseg	= a.current_cseg;
	}

	return *this;
}

int MCMC_Vertex::TotalCodingBits(void)
{
	int nbytes;

	nbytes = 2*sizeof(char)+			// x,y
			sizeof(char)+			//degree
			degree*2*sizeof(char)+
			sizeof(id);

	return nbytes;
}

ofstream & operator>>(MCMC_Vertex &v, ofstream &strm)
{
	(MCMC_MemoryNode)v >> strm;

	// set the space
	strm.write((char*)&(v.Max_Num_Vertex_Degree), sizeof(v.Max_Num_Vertex_Degree));
	// degree
	strm.write((char*)&(v.degree), sizeof(v.degree));
	// coordiniates
	strm.write((char*)&(v.x), sizeof(v.x));
	strm.write((char*)&(v.y), sizeof(v.y));
	for (int i=0; i<v.degree; i++)
		strm.write((char*)&((v.p_contour_segments)[i]),sizeof((v.p_contour_segments)[i]));

	return strm;
}
 
ifstream & operator<<(MCMC_Vertex &s, ifstream &strm)
{
	// set the space
	int temp;

	(MCMC_MemoryNode)s << strm;

	strm.read((char*)&(temp), sizeof(s.Max_Num_Vertex_Degree));
	temp = min(temp, 100);
	
	s.GrowDegreeTo(temp);

	unsigned char degree_temp;
	
	// degree
	s.degree = 0;
	strm.read((char*)&(degree_temp), sizeof(degree_temp));
	if (degree_temp > s.Max_Num_Vertex_Degree)
		s.GrowDegreeTo(s.Max_Num_Vertex_Degree*2);
	s.degree = degree_temp;

	// coordiniates
	strm.read((char*)&(s.x), sizeof(s.x));
	strm.read((char*)&(s.y), sizeof(s.y));
	for (int i=0; i<s.degree; i++)
		strm.read((char*)&((s.p_contour_segments)[i]), sizeof((s.p_contour_segments)[i]));

	return strm;
}

// access
int MCMC_Vertex::GetFirstCSeg(void) // clockwise
{
	assert(degree <= Max_Num_Vertex_Degree);

	if (degree < 1)
		return -1;
	else
	{
		current_cseg = 0;
		return p_contour_segments[current_cseg];
	}
}

int MCMC_Vertex::GetNextCSeg(void)	// clockwise
{
	assert(degree <= Max_Num_Vertex_Degree);
	current_cseg++;
	if (current_cseg >= degree)
		return -1;
	else
		return p_contour_segments[current_cseg];
}

void MCMC_Vertex::GetCSegOfAtomicRegion(mcmcVector<int> &vc_csegs, const int id_region, const int side)	// side 0: right  1:left
{
	int						id_cseg,id_vertex,num=0;
	MCMC_Contour_Segment	*pcseg;
	static int				temp[20];

	vc_csegs.SetRowVector(1);
	vc_csegs.InitValue(-1);

	id_vertex = Id();
	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		if (pcseg->BelongToAtomicRegion(id_region) && pcseg->GetDonotgoRegion()!=id_region)
		{
			int ver1,ver2;
			pcseg->GetVertices(ver1,ver2);
			switch (side)
			{
			case 1: // region at the left side
					if (id_vertex==ver1 && id_region==pcseg->GetLeftAtomicRegion())
						temp[num++] = id_cseg;
					if (id_vertex==ver2 && id_region==pcseg->GetRightAtomicRegion())
						temp[num++] = id_cseg;
			case 0:	// region at the right side
					if (id_vertex==ver1 && id_region==pcseg->GetRightAtomicRegion())
						temp[num++] = id_cseg;
					if (id_vertex==ver2 && id_region==pcseg->GetLeftAtomicRegion())
						temp[num++] = id_cseg;
			}
		}
		id_cseg = GetNextCSeg();
	}
	if (num>0)
	{
		vc_csegs.SetRowVector(num);
		for (int i=0; i<num; i++)
			vc_csegs(i) = temp[i];
	}
}

// side 0: right  1:left
void MCMC_Vertex::GetCSegOfGRegion(mcmcVector<int> &vc_csegs, const int indexid_region, const int side)
{
	int			id_cseg,id_vertex,num=0;
	static int	temp[50];
	MCMC_Contour_Segment	*pcseg;
	bool	bleft_in_gregion,bright_in_gregion;

	vc_csegs.SetRowVector(1);
	vc_csegs.InitValue(-1);

	id_vertex = Id();
	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		// get its left and right atomic regions
		bleft_in_gregion = LayerRegion()->IsInGRegion(pcseg->GetLeftAtomicRegion(), indexid_region);
		bright_in_gregion = LayerRegion()->IsInGRegion(pcseg->GetRightAtomicRegion(), indexid_region);

		if (bleft_in_gregion ^ bright_in_gregion)
		{
			int ver1,ver2;
			pcseg->GetVertices(ver1,ver2);
			switch (side)
			{
			case 1: // region at the left side
					if (id_vertex==ver1 && bleft_in_gregion)
						temp[num++] = id_cseg;
					if (id_vertex==ver2 && bright_in_gregion)
						temp[num++] = id_cseg;
			case 0:	// region at the right side
					if (id_vertex==ver1 && bright_in_gregion)
						temp[num++] = id_cseg;
					if (id_vertex==ver2 && bleft_in_gregion)
						temp[num++] = id_cseg;
			}
		}
		id_cseg = GetNextCSeg();
	}
	if (num>0)
	{
		vc_csegs.SetRowVector(num);
		for (int i=0; i<num; i++)
			vc_csegs(i) = temp[i];
	}
}


// operations
void MCMC_Vertex::Init(void)
{
	Empty();
	// allocated memory
	degree = 0;
	Max_Num_Vertex_Degree = 3;
	p_contour_segments = new int[Max_Num_Vertex_Degree];
	assert(p_contour_segments != NULL);
}

void MCMC_Vertex::Empty(void)
{
	degree = 0;
	// current contour segment
	current_cseg = -1;

	if (p_contour_segments != NULL)
		delete []p_contour_segments;
	p_contour_segments = NULL;
}

void MCMC_Vertex::GrowDegreeTo(const int max_degree)
{
	int *ptemp,i;

	if (Max_Num_Vertex_Degree < max_degree)
	{
		assert(degree <= Max_Num_Vertex_Degree);

		ptemp = new int[max_degree];
		assert(ptemp != NULL);
		for (i=0; i<degree; i++)
			ptemp[i] = p_contour_segments[i];
		delete []p_contour_segments;
		p_contour_segments = ptemp;
		Max_Num_Vertex_Degree = max_degree;
	}
}

int MCMC_Vertex::AddSeg(const int id_cseg)
{
	int i;

	assert(degree <= Max_Num_Vertex_Degree);
	
	if (id_cseg<0)
		return 0;	// valid id?
	for (i=0; i<degree; i++)
		if (p_contour_segments[i]==id_cseg)	// segment exists?
			return 0;
	if (degree>=Max_Num_Vertex_Degree)
		GrowDegreeTo(Max_Num_Vertex_Degree*2);

	assert(degree < Max_Num_Vertex_Degree);	// add
	assert (p_contour_segments!=NULL);
	if (p_contour_segments!=NULL)
	{
		p_contour_segments[degree] = id_cseg;
		degree++;
	}
	return 1;
}

int MCMC_Vertex::DeleteCSeg(const int id_cseg)	// delete segment index from this vertex
{
	int i,j;
	i = GetFirstCSeg();		// get the first contour segment
	while (i >= 0)
	{
		if (i == id_cseg)
		{
			assert(degree <= Max_Num_Vertex_Degree);

			for (j=current_cseg; j<degree-1; j++)
				p_contour_segments[j] = p_contour_segments[j+1];
			degree--;		// degree decrease
			if (degree == 0)// become dangling vertex
				LayerRegion()->DeleteVertexFromArray(Id());

			return 1;		// suceed
		}
		i = GetNextCSeg();	// get the next contour segment
	}
	return 0;				// fail
}

int	MCMC_Vertex::GetNode(void)
{
	if (degree <= 0)
		return -1;
	else
	{
		MCMC_Contour_Segment *pcseg;

		assert(degree <= Max_Num_Vertex_Degree);

		pcseg = LayerRegion()->GetCSeg(p_contour_segments[0]);
		assert(pcseg != NULL);
		if (pcseg->end[0] == Id())
			return pcseg->node[0];
		else
		{
			assert(pcseg->end[1] == Id());
			return pcseg->node[1];
		}
	}
}

// total degree with self-closed contour counted as two
int MCMC_Vertex::DegreeExt(void)
{
	int id_cseg,count=0;
	MCMC_Contour_Segment *pcseg;

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		if (pcseg->SelfClosed())
			count += 2;
		else
			count ++;

		id_cseg = GetNextCSeg();
	}

	return count;
}


//////////////////////////////////////////////////////////////
// MCMC Countour Segment
/////////////////////////////////////////////////////////////
// construction & destruction
MCMC_Contour_Segment::MCMC_Contour_Segment()
{
	player_region = NULL;

	Empty();
}

MCMC_Contour_Segment::~MCMC_Contour_Segment()
{
	player_region = NULL;
	Empty();
}

void MCMC_Contour_Segment::Empty(void)
{
	boundry_of[0] = -1;	boundry_of[1] = -1;
	//	two extreme vertices
	end[0] = -1;	end[1] = -1;
	//	two extreme nodes
	node[0] = -1;	node[1] = -1;
	extn[0] = -1;	extn[1] = -1;
	region_donotgo = -1;			// the region that should not go along this segment
}

// operators
ofstream & operator>>(MCMC_Contour_Segment &x, ofstream &strm)
{
	(MCMC_MemoryNode)x >> strm;

	for (int i=0; i<2; i++)
	{
		strm.write((char*)&((x.boundry_of)[i]), sizeof((x.boundry_of)[i]));
		strm.write((char*)&((x.end)[i]), sizeof((x.end)[i]));
		strm.write((char*)&((x.node)[i]), sizeof((x.node)[i]));
		strm.write((char*)&((x.extn)[i]), sizeof((x.extn)[i]));
	}
	strm.write((char*)&(x.region_donotgo), sizeof(x.region_donotgo));

	return strm;
}

ifstream &	operator<<(MCMC_Contour_Segment &x, ifstream &strm)
{
	(MCMC_MemoryNode)x << strm;

	for (int i=0; i<2; i++)
	{
		strm.read((char*)&((x.boundry_of)[i]), sizeof((x.boundry_of)[i]));
		strm.read((char*)&((x.end)[i]), sizeof((x.end)[i]));
		strm.read((char*)&((x.node)[i]), sizeof((x.node)[i]));
		strm.read((char*)&((x.extn)[i]), sizeof((x.extn)[i]));
	}
	strm.read((char*)&(x.region_donotgo), sizeof(x.region_donotgo));

	return strm;
}

MCMC_Contour_Segment & MCMC_Contour_Segment::operator=(const MCMC_Contour_Segment &a)
{
	if (this != &a)
	{
		boundry_of[0]	= a.boundry_of[0];
		boundry_of[1]	= a.boundry_of[1];
		end[0]		= a.end[0];
		end[1]		= a.end[1];
		node[0]		= a.node[0];
		node[1]		= a.node[1];
		extn[0]		= a.extn[0];
		extn[1]		= a.extn[1];
		direction	= a.direction;
		current_node= a.current_node;
	}

	return *this;
}

int MCMC_Contour_Segment::TotalCodingBits(void)
{
	int		nbytes;
	int			id_node;
	MCMC_Node	*pnode;
	double		dbytes=0.0;

	nbytes = 3*sizeof(char)+	//labels
		    3*sizeof(char)+		// ends
			3*sizeof(char);		// nodes

	id_node = GetFirstNode();
	while (id_node >= 0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		dbytes += pnode->TotalCodingBits();

		id_node = GetNextNode();
	}

	nbytes += (int)dbytes;
	return nbytes;
}


// access
void MCMC_Contour_Segment::GetVertices(int &id_ver1, int &id_ver2)
{
	id_ver1 = end[0];
	id_ver2 = end[1];
}

int MCMC_Contour_Segment::GetTheOtherVertex(const int id_ver1)
{
	if (id_ver1 == end[0])
		return end[1];
	else
	{
		if (id_ver1 == end[1])
			return end[0];
		else
			return -1;
	}
}

bool MCMC_Contour_Segment::BelongToAtomicRegion(const int id_region)
{
	if (boundry_of[0]==id_region || boundry_of[1]==id_region)
		return true;
	else
		return false;
}

bool MCMC_Contour_Segment::IsOfAtomicRegions(const int id_region1, const int id_region2)
{
	int i1 = GetLeftAtomicRegion();
	int i2 = GetRightAtomicRegion();
	if ((id_region1==i1 && id_region2==i2) || (id_region1==i2 && id_region2==i1))
		return true;
	else
		return false;
}

// operations
int MCMC_Contour_Segment::ReplaceAtomicRegionId(const int id_rgnS, const int id_rgnD)
{
	if (boundry_of[0] == id_rgnS)
	{
		boundry_of[0] = id_rgnD;
		return 1;
	}
	else
	{
		if (boundry_of[1] == id_rgnS)
		{
			boundry_of[1] = id_rgnD;
			return 1;
		}
	}
	return 0;
}

void MCMC_Contour_Segment::SwapVertices(void)
{
	MCMC_Node	*pnode;
	int			id_node;

	id_node = GetFirstNode();
	while (id_node >= 0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);

		id_node = GetNextNode();
		SwapInt(pnode->link[0],pnode->link[1]);
	}
	// switch region labels
	SwapInt(boundry_of[0], boundry_of[1]);
	// switch extreme vertices
	SwapInt(end[0], end[1]);
	// switch two extreme nodes
	SwapInt(node[0],node[1]);
	// switch two extentions
	SwapInt(extn[0],extn[1]);
}

int MCMC_Contour_Segment::GetFirstNode(void)
{
	current_node = node[0];
	return current_node;
}

int MCMC_Contour_Segment::GetNextNode(void)
{
	if (node[0]==node[1] && extn[0]==extn[1] &&
		node[0]==extn[0])
		return -1;

	if (current_node<0)
		return -1;

	if (current_node == node[0])
	{
		// current pointer is with node[0];
		current_node = extn[0];
		return current_node;
	}

	if (current_node == extn[1])
	{
		if (node[0] != node[1])
		{
			current_node = node[1];
			return current_node;
		}
		else
			return -1;
	}
	if (current_node == node[1])
	{
		current_node = -1;
		return current_node;
	}
	MCMC_Node *pnode;
	pnode = LayerRegion()->GetNode(current_node);
	if (pnode == NULL)
		return -1;
	// ?? forget why?
	if (pnode->link[0] == extn[0])
		current_node = pnode->link[1];
	else
		current_node = pnode->link[0];
	return current_node;
}

// in case this segment is a self-closed segment
int	MCMC_Contour_Segment::GetNextNodeExtn(void)
{
	if (!SelfClosed())
		return GetNextNode();
	else
	{
		if (current_node==-2)
			current_node = -1;
		else if (current_node==extn[1])
			current_node = -2;
		else
			GetNextNode();

		if (current_node==-2)
			return node[1];
		else
			return current_node;
	}
}

int MCMC_Contour_Segment::GetLastNode(void)
{
	current_node = node[1];
	return current_node;
}

int MCMC_Contour_Segment::GetPrevNode(void)
{
	if (current_node<0)
		return -1;

	if (current_node == node[1])
	{
		// current pointer is with node[1];
		current_node = extn[1];
		return current_node;
	}

	if (current_node == extn[0])
	{
		if (node[1] != node[0])
		{
			current_node = node[0];
			return current_node;
		}
		else
			return -1;
	}

	if (current_node == node[0])
	{
		current_node = -1;
		return current_node;
	}

	MCMC_Node *pnode;
	pnode = LayerRegion()->GetNode(current_node);
	assert(pnode != NULL);
	if (pnode->link[1] == extn[1])
		current_node = pnode->link[0];
	else
		current_node = pnode->link[1];
	return current_node;
}

// draw this segment to image
bool MCMC_Contour_Segment::PointOfLine(const int x, const int y)
{
	int			id_node;
	MCMC_Node	*pnode;

	id_node = GetFirstNode();
	while (id_node >= 0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		if (pnode->GetY()==x && pnode->GetX()==y)
			return true;					// found one

		id_node = GetNextNode();
	}
	return false;
}

// get the length of this segment
double MCMC_Contour_Segment::GetLength(void)
{
	int			id_node;
	double		dLength=0.0,x1=0.0,y1=0.0,x2,y2;
	MCMC_Node	*pnode;
	bool		bFirst=true;	

	id_node = GetFirstNode();
	while (id_node>=0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		if (bFirst)
			bFirst = false;
		else
		{
			x2 = pnode->GetY();
			y2 = pnode->GetX();
			dLength += Distance(x1,y1,x2,y2);
		}
		x1 = pnode->GetY();
		y1 = pnode->GetX();

		id_node = GetNextNode();
	}
	
	return dLength;
}

int	MCMC_Contour_Segment::GetTheOtherRegionId(const int id_region1)
{
	if (id_region1 == GetLeftAtomicRegion())
		return GetRightAtomicRegion();
	else
	{
		if (id_region1!=GetRightAtomicRegion())
		{
			cout << "GetTheOtherRegionId" << endl;
			XPT_OutOfBound *pxpt;
			pxpt = new XPT_OutOfBound("Out of nodes");

			throw pxpt;
		}
		return GetLeftAtomicRegion();
	}
}

int	MCMC_Contour_Segment::GetNumNode(void)
{
	int id_node,num=0;

	id_node = GetFirstNode();
	while (id_node>=0)
	{
		num++;
		id_node = GetNextNode();
	}
	return num;
}

bool MCMC_Contour_Segment::ShowOnPixel(const int id_region, const int id_node)
{
	if (id_node==node[0] || id_node==node[1] ||
		id_node==extn[0] || id_node==extn[1])
		return true;

	int id_node_temp;
	MCMC_AtomicRegion *pregion;
	MCMC_Node *pnode;

	pregion = LayerRegion()->GetAtomicRegion(id_region);
	pnode = LayerRegion()->GetNode(id_node);
	if (pregion->Direction())
		id_node_temp = id_node;
	else
		id_node_temp = pnode->link[1];

	pnode = LayerRegion()->GetNode(id_node_temp);
	if (pnode != NULL)
	{
		if (GetLeftAtomicRegion()==id_region)
			return true;
		else if (GetRightAtomicRegion()==id_region)
			return true;
		else
			return false;
	}
	else
		return true;
}

void MCMC_Contour_Segment::DeleteNodesInBetween(const int id_node1, const int id_node2)
{
	int			id_node,id_temp;
	MCMC_Node	*pnode1,*pnode2,*pnode;
	int			nleft_node=0,nright_node=0;
	McPoint2D	pt2d,ptLnBg,ptLnEd;

	pnode1 = LayerRegion()->GetNode(id_node1);
	pnode2 = LayerRegion()->GetNode(id_node2);
	if (pnode1==NULL || pnode2==NULL || id_node1==id_node2)
		return;
	// change coordinates system
	ptLnBg.Set(pnode1->GetX(),pnode1->GetY());
	ptLnEd.Set(pnode2->GetX(),pnode2->GetY());

	// just delete them
	RestoreCurrent(id_node1);
	id_node = GetNextNodeExtn();
	while (id_node != id_node2)
	{
		id_temp = id_node;
		pnode = LayerRegion()->GetNode(id_temp);
		pt2d.Set(pnode->GetX(),pnode->GetY());
		if (PointOnLeftSide(pt2d,ptLnBg,ptLnEd))
			nleft_node++;
		else
			nright_node++;

		id_node = GetNextNodeExtn();
		LayerRegion()->DeleteNodeFromArray(id_temp);
	}
	pnode1->link[0] = id_node2;
	pnode2->link[1] = id_node1;

	if (id_node1==node[0])
		extn[0] = id_node2;
	if (id_node2==node[1])
		extn[1] = id_node1;
}

void MCMC_Contour_Segment::InsertNodesInBetween(const int id_node1, const int id_node2, const bool bjust_one)
{
	MCMC_Node	*pnode1,*pnode2,*pnode_new;
	MCMC_Node	*pnode_temp1,*pnode_temp2;
	int			id_node_new,num_new_node;
	int			id_temp1,id_temp2,i;
	double		dist;

	pnode1 = LayerRegion()->GetNode(id_node1);
	pnode2 = LayerRegion()->GetNode(id_node2);
	if (pnode1==NULL || pnode2==NULL || id_node1==id_node2)
		return;

	dist = Distance2(pnode1->GetX(),pnode1->GetY(),pnode2->GetX(),pnode2->GetY());
	if (dist<MINI_DIST1*1.5)
		return;

	if (bjust_one)
		num_new_node = 1;
	else
		num_new_node = (int)round(dist/MINI_DIST1)-1;

	id_temp1 = id_node1;
	id_temp2 = id_node2;
	for (i=0; i<num_new_node; i++)
	{
		pnode_temp1 = LayerRegion()->GetNode(id_temp1);
		pnode_temp2 = LayerRegion()->GetNode(id_temp2);
		id_node_new = LayerRegion()->NewNode();
		pnode_new = LayerRegion()->GetNode(id_node_new);
		assert(pnode_new != NULL);
		pnode_new->link[1] = id_temp1;
		pnode_new->link[0] = id_temp2;
		pnode_temp1->link[0] = id_node_new;
		pnode_temp2->link[1] = id_node_new;

		dist = Distance2(pnode_temp1->GetX(),pnode_temp1->GetY(),pnode_temp2->GetX(),pnode_temp2->GetY());
		if (dist>1.9*MINI_DIST1 && !bjust_one)
		{
			pnode_new->SetX(pnode_temp1->GetX()+(pnode_temp2->GetX()-pnode_temp1->GetX())/dist*MINI_DIST1);
			pnode_new->SetY(pnode_temp1->GetY()+(pnode_temp2->GetY()-pnode_temp1->GetY())/dist*MINI_DIST1);
		}
		else
		{
			pnode_new->SetX(pnode_temp1->GetX()/2.0+pnode_temp2->GetX()/2.0);
			pnode_new->SetY(pnode_temp1->GetY()/2.0+pnode_temp2->GetY()/2.0);
		}
		id_temp1 = id_node_new;
	}
	if (num_new_node>0)
	{
		if (id_node1==node[0])
			extn[0] = pnode1->link[0];
		if (id_node2==node[1])
			extn[1] = pnode2->link[1];
	}
}


void MCMC_Contour_Segment::Simplify(const int id_node1, const int id_node2)
{
	int			id_node;
	double		dLen;
	MCMC_Node	*pnode,*pnode1,*pnode2;

	pnode1 = LayerRegion()->GetNode(id_node1);
	pnode2 = LayerRegion()->GetNode(id_node2);
	if (pnode1==NULL || pnode2==NULL || id_node1==id_node2)
		return;

	dLen = 0.0;
	RestoreCurrent(id_node1);
	pnode2 = pnode1;
	id_node = GetNextNode();
	while (id_node>=0 && id_node!=id_node2)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		dLen += Distance(pnode2->GetY(),pnode2->GetX(),pnode->GetY(),pnode->GetX());

		if (dLen > MINI_DIST1)
			break;
		id_node = GetNextNode();
		pnode2 = pnode;
	}
	if (dLen > MINI_DIST1)
	{
		DeleteNodesInBetween(id_node1, id_node);
		Simplify(id_node,id_node2);
	}
	else
		DeleteNodesInBetween(id_node1, id_node2);

	// insert some nodes here if the distance between two nodes are too large
	int			id_node_new,id1,id2;
	MCMC_Node	*pnode_new;

	RestoreCurrent(id_node1);
	id1 = id_node1;
	id2 = GetNextNode();
	while (id2>=0 && id2!=id_node2)
	{
		pnode1 = LayerRegion()->GetNode(id1);
		pnode2 = LayerRegion()->GetNode(id2);
		if (Distance(pnode2->GetY(),pnode2->GetX(),pnode1->GetY(),pnode1->GetX())>MINI_DIST1/2.0*3.0)
		{
			id_node_new = LayerRegion()->NewNode();
			pnode_new = LayerRegion()->GetNode(id_node_new);
			assert(pnode_new != NULL);

			pnode_new->SetX((pnode2->GetX()+pnode1->GetX())/2.0);
			pnode_new->SetY((pnode2->GetY()+pnode1->GetY())/2.0);
			pnode_new->link[0] = id2;
			pnode_new->link[1] = id1;
			pnode1->link[0] = id_node_new;
			pnode2->link[1] = id_node_new;

			if (id1==node[0])
				extn[0] = id_node_new;
			if (id2==node[1])
				extn[1] = id_node_new;
		}

		id1 = id2;
		id2 = GetNextNode();
	}
}

void MCMC_Contour_Segment::Simplification(const int id_node1, const int id_node2, mcmcMatrix<double> &distances)
{
	int			id_cur_temp,i,id_node_min,id_node;
	MCMC_Node	*pnode1,*pnode2,*pnode;
	McPoint2D	pt2d,ptLnBg,ptLnEd;
	double		dMinDis;

	if (id_node1<0 || id_node2<0)
		return;

	// set initial values
	id_node_min = -1;
	dMinDis = MINI_SIMP_DIST;

	// save current node
	id_cur_temp = KeepCurrent();
	RestoreCurrent(id_node1);
	pnode1 = LayerRegion()->GetNode(id_node1);
	assert(pnode1 != NULL);
	ptLnBg.Set(pnode1->GetY(),pnode1->GetX());

	pnode2 = LayerRegion()->GetNode(id_node2);
	assert(pnode2 != NULL);
	ptLnEd.Set(pnode2->GetY(),pnode2->GetX());

	i = 0;
	id_node = GetNextNode();
	while (id_node>=0 && id_node!=id_node2)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		pt2d.Set(pnode->GetY(),pnode->GetX());
		if (id_node1 != id_node2)
			distances(i) = PerDistance2(pt2d,ptLnBg,ptLnEd);	// common nodes
		else
			distances(i) = Distance(pt2d,ptLnBg);			// the vertex has degree two

		if (distances(i)>dMinDis)
		{
			dMinDis = distances(i);
			id_node_min = id_node;
		}
		// keep the corner node anyway
		if (LayerRegion()->OnImageCorner(pt2d.x,pt2d.y))
		{
			dMinDis = 1.0;
			id_node_min = id_node;
			break;
		}
		i++;
		id_node = GetNextNode();
	}
	if (id_node_min>=0)
	{
		Simplification(id_node1,id_node_min,distances);
		Simplification(id_node_min,id_node2,distances);
	}
	else
		Simplify(id_node1,id_node2);

	// set back
	RestoreCurrent(id_cur_temp);
}

void MCMC_Contour_Segment::DeleteVeryColseNodes(void)
{
	int			id_node1,id_node,id_node2,num;
	MCMC_Node	*pnode1,*pnode,*pnode2;

	num = GetNumNode();

	id_node1 = GetFirstNode();
	id_node = GetNextNodeExtn();
	id_node2 = GetNextNodeExtn();
	while (id_node1>=0 && id_node>=0 && id_node2>=0)
	{
		int current_temp = current_node;

		pnode1 = LayerRegion()->GetNode(id_node1);
		pnode = LayerRegion()->GetNode(id_node);
		pnode2 = LayerRegion()->GetNode(id_node2);

		bool bDel=false;
		if (Distance(pnode1->GetX(),pnode1->GetY(),pnode->GetX(),pnode->GetY())<2.1
			|| Distance(pnode2->GetX(),pnode2->GetY(),pnode->GetX(),pnode->GetY())<2.1)
			bDel = true;
		// don't delete any node on the corners
		if (LayerRegion()->OnImageCorner(pnode->GetX(),pnode->GetY()))
			bDel = false;

		if (bDel)
		{
			DeleteNodesInBetween(id_node1, id_node2);
			current_node = current_temp;
			// keep id_node1
			id_node = id_node2;
			id_node2 = GetNextNodeExtn();
		}
		else
		{
			id_node1 = id_node;
			id_node = id_node2;
			id_node2 = GetNextNodeExtn();
		}
	}
}

void MCMC_Contour_Segment::Smooth(void)
{
	int				id_node1,id_node2;
	mcmcMatrix<double>	distances;
	bool			bShouldSmooth=false;

	// insert nodes if the distance between any two nodes is too large
	id_node1 = GetFirstNode();
	id_node2 = GetNextNode();
	while (id_node1>=0 && id_node2>=0)
	{
		InsertNodesInBetween(id_node1,id_node2);

		id_node1 = id_node2;
		id_node2 = GetNextNode();
	}

	// delete nodes if the distance between any two nodes is too small
	distances.SetDimension(GetNumNode(),1);
	Simplification(node[0], node[1],distances);

	DeleteVeryColseNodes();

	// insert nodes if the distance between any two nodes is too large
	id_node1 = GetFirstNode();
	id_node2 = GetNextNodeExtn();
	while (id_node1>=0 && id_node2>=0)
	{
		InsertNodesInBetween(id_node1,id_node2,false);

		id_node1 = id_node2;
		id_node2 = GetNextNodeExtn();
	}
}


bool EqualValue(const double a1, const double a2)
{
	if (fabs(a1-a2)<1e-5)
		return true;
	else
		return false;
}



int	MCMC_Contour_Segment::GetNextNodeToVertex(const int id_ver)
{
	if (id_ver==end[0])
		return extn[0];
	else
	{
		assert(id_ver==end[1]);
		return extn[1];
	}
}

void MCMC_Contour_Segment::Draw(Raster<int> &image,const int value)
{
	int			node1,node2,x1,y1,x2,y2,i,i1,j;
	MCMC_Node	*pnode1,*pnode2;
	bool		bcircle_back;

	// self-closed contour segment and come back to the start point?
	bcircle_back = false;

	if (!Valid())
		return;

	node1 = GetFirstNode();
	node2 = GetNextNode();
	while (node1>=0 && node2>=0)
	{
		// get coordinates
		pnode1 = LayerRegion()->GetNode(node1);
		assert(pnode1 != NULL);
		pnode2 = LayerRegion()->GetNode(node2);
		assert(pnode2 != NULL);
		x1 = IntClose(pnode1->GetY()); y1 = IntClose(pnode1->GetX());
		x2 = IntClose(pnode2->GetY()); y2 = IntClose(pnode2->GetX());

		double theta;
		if (x1!=x2 && abs(x2-x1)>abs(y2-y1))
		{
			theta = atan2(double(y2-y1),double(x2-x1));
			if (x2<x1)	i1 = -1;
			else		i1 = 1;
			for (i=x1; i!=x2; i=i+i1)
			{
				j = y1+IntClose(tan(theta)*(i-x1));
				if (j>=0 && j<=LayerRegion()->height && i>=0 && i<=LayerRegion()->width)
				{
					if (image.Valid(j,i))
						image(j,i) = value;
				}
			}
		}
		else
		{
			theta = atan2(double(y2-y1),double(x2-x1));
			if (y2<y1)	i1 = -1;
			else		i1 = 1;
			for (j=y1; j!=y2; j=j+i1)
			{
				i = x1+IntClose(1.0/tan(theta)*(j-y1));
				if (j>=0 && j<=LayerRegion()->height && i>=0 && i<=LayerRegion()->width)
				{
					if (image.Valid(j,i))
						image(j,i) = value;
				}
			}
		}
		// self-closed contour segment and come back to the start point?
		if (bcircle_back)
			break;

		node1 = node2;
		node2 = GetNextNode();

		if (SelfClosed() && node2<0 && node1>=0)
		{
			node2 = GetFirstNode();
			bcircle_back = true;
		}
	}
}

bool MCMC_Contour_Segment::IsBoundaryOfGRegion(const int indexid)
{
	bool	bleft,bright;

	bleft = LayerRegion()->IsInGRegion(GetLeftAtomicRegion(), indexid);
	bright = LayerRegion()->IsInGRegion(GetRightAtomicRegion(), indexid);

	return (bleft ^ bright);
}

bool MCMC_Contour_Segment::OnlyTwoNodes(const int id_node1, const int id_node2)
{
	if ((node[0]==id_node1 && node[1]==id_node2 && extn[0]==node[1]) ||
		(node[1]==id_node1 && node[0]==id_node2 && extn[0]==node[1]))
		return true;
	else
		return false;
}

void MCMC_Contour_Segment::ChangeStartOfRelatedRegion()
{
	MCMC_AtomicRegion *pregion;

	pregion = LayerRegion()->GetAtomicRegion(GetLeftAtomicRegion());
	if (pregion!=NULL && pregion->GetStartCSeg()==Id())
		pregion->UseDifferentStartCSeg();

	pregion = LayerRegion()->GetAtomicRegion(GetRightAtomicRegion());
	if (pregion!=NULL && pregion->GetStartCSeg()==Id())
		pregion->UseDifferentStartCSeg();
}

void MCMC_Contour_Segment::UpdateTopologyRecord(void)
{
	MCMC_Vertex	*pver;
	int			id_node1,id_node2;
	MCMC_Node	*pnode1,*pnode2;

	pver = LayerRegion()->GetVertex(end[0]);
	assert(pver != NULL);
	LayerRegion()->mx_label_id_vertices(pver->GridRow(), pver->GridCol()) = pver->Id();

	id_node1 = GetFirstNode();
	id_node2 = GetNextNode();

	while (id_node1 >= 0)
	{
		pnode1 = LayerRegion()->GetNode(id_node1);
		pnode2 = LayerRegion()->GetNode(id_node2);

		// the node position
		LayerRegion()->mx_label_id_nodes(pnode1->GridRow(), pnode1->GridCol()) = id_node1;
		if (pnode1!=NULL && pnode2!=NULL)
			LayerRegion()->cseg_on_grid.CSegId(pnode1->GridRow(), pnode1->GridCol(), pnode2->GridRow(), pnode2->GridCol()) = Id();

		id_node1 = id_node2;
		id_node2 = GetNextNode();
	}
	if (node[0] == node[1])
	{
		// self-closed contour
		pnode1 = LayerRegion()->GetNode(node[1]);
		pnode2 = LayerRegion()->GetNode(extn[1]);
		assert(pnode1!=NULL && pnode2!=NULL);
		LayerRegion()->cseg_on_grid.CSegId(pnode1->GridRow(), pnode1->GridCol(), pnode2->GridRow(), pnode2->GridCol()) = Id();
	}
	else
	{
		pver = LayerRegion()->GetVertex(end[1]);
		assert(pver != NULL);
		LayerRegion()->mx_label_id_vertices(pver->GridRow(), pver->GridCol()) = pver->Id();
	}
}


//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
	// construction & destruction
MCMC_StartingSeg::MCMC_StartingSeg()
{
	pregion1 = NULL;
	pregion2 = NULL;
	num_seg = 0;
}

MCMC_StartingSeg::~MCMC_StartingSeg()
{
	pregion1 = NULL;
	pregion2 = NULL;
}

// operations
void MCMC_StartingSeg::Set(MCMC_AtomicRegion &region1, MCMC_AtomicRegion &region2, const bool bAddSegs)
{
	int						id_cseg,i;
	MCMC_Contour_Segment	*pcseg;
	mcmcVector<int>				vc_cseg;

	pregion1 = &region1;
	pregion2 = &region2;

	// initilization
	sign.SetDimension(region1.LayerRegion()->MaxCSegSpace(),1);
	sign.InitValue(0);
	num_seg = 0;
	id_start_temp = -1;

	if (bAddSegs)
	{
		// store every segments that belong to region1 or region2
		region1.LayerRegion()->AllCSegToMatrix(vc_cseg);
		for (i=0; i<vc_cseg.size(); i++)
		{
			id_cseg = vc_cseg(i);
			pcseg = region1.LayerRegion()->GetCSeg(id_cseg);
			assert(pcseg!=NULL);
			if (pcseg->GetLeftAtomicRegion()==region1.Id() || pcseg->GetLeftAtomicRegion()==region2.Id() ||
				pcseg->GetRightAtomicRegion()==region1.Id() || pcseg->GetRightAtomicRegion()==region2.Id())
			{
				sign(id_cseg) = 1;
				num_seg++;
			}
		}
	}
}

void MCMC_StartingSeg::Delete(const int id_cseg)
{
	assert(sign(id_cseg)!=0);
	sign(id_cseg) = 0;
	num_seg--;
}

void MCMC_StartingSeg::Keep(const int id_cseg)
{
	id_start_temp = id_cseg;
}

void MCMC_StartingSeg::Add(const int id_cseg)
{
	if (sign.size()==0)
		return;
	if (sign(id_cseg)==0)
	{
		sign(id_cseg) = 1;
		num_seg++;
	}
	id_start_temp = id_cseg;
}

int MCMC_StartingSeg::GetStartingSeg(void)
{
	int	id_cseg,num;
	int id_temp = pregion1->GetStartCSeg();
	int MAX_TRY,ntry=0;
	LinkList<int>	list_cseg;

	MAX_TRY = pregion1->LayerRegion()->MaxCSegSpace();
	// region 2 is a hole of region 1
	if (pregion1->IsHole(pregion2->Id()))
		return pregion1->GetStartCSeg();

	if (id_start_temp<0)
		id_start_temp = pregion1->GetStartCSeg();

	pregion1->ChangeStartCSeg(id_start_temp);
	num = 0;
	id_cseg = pregion1->GetFirstCSeg();
	while (id_cseg>=0 && num<MAX_TRY)
	{
		num++;
		id_cseg = pregion1->GetNextCSeg();
	}
	if (num == num_seg)
	{
		pregion1->ChangeStartCSeg(id_temp);
		return id_start_temp;
	}

	MCMC_Boundary	*pBoundary;
	int				*start_temp;

	pBoundary = new MCMC_Boundary[num_seg];
	start_temp = new int[num_seg];

	assert(pBoundary != NULL);

	int i=0,j=0;
	while (num_seg>0)
	{
		for (j=0; j<sign.size(); j++)
		{
			if (sign(j)!=0)
				break;
		}
		pregion1->ChangeStartCSeg(j);
		start_temp[i] = pregion1->GetStartCSeg();

		id_cseg = pregion1->GetFirstCSeg();
		ntry = 0;
		list_cseg.Empty();
		while (id_cseg>=0 && ntry<MAX_TRY)
		{
			ntry++;
			list_cseg.Insert(id_cseg);
			id_cseg = pregion1->GetNextCSeg();
		}
		if (ntry<MAX_TRY)
		{
			// the starting segment is valid
			bool bGo;
			bGo = list_cseg.GetFirst(id_cseg);
			while (bGo)
			{
				if (sign(id_cseg)!=0)
					num_seg--;
				sign(id_cseg) = 0;
				bGo = list_cseg.GetNext(id_cseg);
			}
			pBoundary[i].Set(*pregion1);
		}
		else
		{
			num_seg -= 1;
			sign(pregion1->GetStartCSeg()) = 0;
		}
		i++;
	}

	double darea=0;
	int	   new_start;

	for (j=0; j<i; j++)
	{
		if (pBoundary[j].Area()>darea)
		{
			new_start = start_temp[j];
			darea = pBoundary[j].Area();
		}
	}

	delete []start_temp;
	delete []pBoundary;

	return new_start;
}

//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
// region statististics

// construction & destruction
MCMC_Region_Stat::MCMC_Region_Stat()
{
	pmodel = NULL;
}

MCMC_Region_Stat::~MCMC_Region_Stat()
{
	DeleteModel();
}

void MCMC_Region_Stat::DeleteModel(void)
{
	if (pmodel!=NULL)
		delete pmodel;
	pmodel = NULL;
}

// operators
ofstream & operator>>(MCMC_Region_Stat &x, ofstream &strm)
{
	if (x.pmodel!=NULL)
		(x.pmodel)->SaveToFile(strm);
	else
	{
		int temp=-1;
		strm.write((char*)&temp, sizeof(temp));
	}

	return strm;
}

ifstream & operator<<(MCMC_Region_Stat &x, ifstream &strm)
{
	int type;

	x.DeleteModel();

	strm.read((char*)&type, sizeof(type));
	if (type<0)
		x.pmodel = NULL;
	else
	{
		x.pmodel = CreateModel(type);
		(x.pmodel)->ReadFromFile(strm);
	}

	return strm;
}

// operation
MCMC_Region_Stat & MCMC_Region_Stat::operator=(MCMC_Region_Stat &x)
{
	if (this != &x && x.Model()!=NULL)
	{
		// free the previous model first
		DeleteModel();
		// create a same model as pmodel1
		pmodel = CreateModel(x.pmodel->Type());
		CopyModel(pmodel,x.pmodel);
	}

	return *this;
}

void MCMC_Region_Stat::SetStat(MCMC_Model *pmodel1)
{
	if (pmodel1!=NULL && pmodel1!=Model())
	{
		// free the previous model first
		DeleteModel();
		pmodel = pmodel1->CopyTo();
	}
}

double MCMC_Region_Stat::GetLogPostProb(const int row_pixel, const int col_pixel)
{
	if (pmodel==NULL)
		return 0.0;
	else
		return pmodel->GetLogPostProb(row_pixel, col_pixel);
}

//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
// sorted node in a region

// construction & destruction
MCMC_NodeVertex::MCMC_NodeVertex()
{
	id_node = -1;
	id_ver = -1;
	id_cseg = -1;
}

MCMC_NodeVertex::~MCMC_NodeVertex()
{
	id_node = -1;
	id_ver = -1;
	id_cseg = -1;
}

// operators
MCMC_NodeVertex& MCMC_NodeVertex::operator=(const MCMC_NodeVertex& x)
{
	if (&x != this)
	{
		id_node = x.id_node;
		id_ver  = x.id_ver;
		id_cseg = x.id_cseg;
		ofOverall = x.ofOverall;
		ofStat  = x.ofStat;
		ofSmooth = x.ofSmooth;
	}
	return *this;
}

bool MCMC_NodeVertex::operator==(const MCMC_NodeVertex& x)
{
	if (id_node == x.id_node)
	{
		assert(id_ver == x.id_ver);
		return true;
	}
	else
		return false;
}

//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////

// construction & destruction
MCMC_SortedNodes::MCMC_SortedNodes()
{
	num_node = 0;
	psorted_nodes = NULL;
	ppcontour_nodes_holes = NULL;
}

MCMC_SortedNodes::~MCMC_SortedNodes()
{
	Empty();
}

void MCMC_SortedNodes::Empty(void)
{
	num_node = 0;

	if (ppcontour_nodes_holes != NULL)
		delete []ppcontour_nodes_holes;
	ppcontour_nodes_holes = NULL;

	if (psorted_nodes != NULL)
		delete []psorted_nodes;
	psorted_nodes = NULL;
}

// operator
MCMC_SortedNodes & MCMC_SortedNodes::operator=(MCMC_SortedNodes &a)
{
	if (this != &a)
	{
		int i,id_region_hole;

		dArea = a.dArea;
		num_node = a.num_node;
		index_move = a.index_move;
		if (psorted_nodes != NULL)
			delete []psorted_nodes;
		psorted_nodes = NULL;
		if (num_node>0)
		{
			psorted_nodes = new MCMC_NodeVertex[num_node];
			for (i=0; i<num_node; i++)
				psorted_nodes[i] = a.psorted_nodes[i];
		}
		if (ppcontour_nodes_holes != NULL)
			delete []ppcontour_nodes_holes;
		ppcontour_nodes_holes = NULL;
		if (NumHoles()>0)
		{
			ppcontour_nodes_holes = new MCMC_PSortedNodes[NumHoles()];
			// add the boundary of its hole
			i = 0;
			for (int j=0; j<NumHoles(); j++)
			{
				id_region_hole = GetHole(j);

				ppcontour_nodes_holes[i++] = LayerRegion()->GetSortedNodes(id_region_hole);
			}
		}
	}

	return *this;
}

ofstream & operator>>(MCMC_SortedNodes &a, ofstream &strm)
{
	strm.write((char*)&(a.dArea),sizeof(a.dArea));

	return strm;
}

ifstream & operator<<(MCMC_SortedNodes &a, ifstream &strm)
{
	strm.read((char*)&(a.dArea),sizeof(a.dArea));

	return strm;
}

double MCMC_SortedNodes::ComputeArea(void)
{
	double		*px,*py;
	MCMC_Node	*pnode;
	int			i,id_node;

	dArea = 0.0;
	num_node = NumNodes();
	if (num_node<=0)
		return Area();

	px = new double[num_node];
	py = new double[num_node];
	i = 0;
	id_node = GetFirstNode();
	while (id_node>=0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		assert(pnode != NULL);
		px[i] = pnode->GetX();
		py[i++] = pnode->GetY();

		id_node = GetNextNode();
	}
	// compute the total area
	dArea = area(px,py,num_node);

	delete []px;
	delete []py;

	return dArea;
}

// gathering all the nodes in clockwise order
void MCMC_SortedNodes::SortAllNodes(void)
{
	bool		bClockwise;
	int			i,id_node;

	num_node = NumNodes();
	if (psorted_nodes != NULL)
		delete []psorted_nodes;
	psorted_nodes = new MCMC_NodeVertex[num_node];

	assert(psorted_nodes != NULL);
	// node order in clockwise or counterclockwise
	SetDirection();
	bClockwise = Clockwise();

	i = 0;
	id_node = GetFirstNode();
	while (id_node>=0)
	{
		psorted_nodes[i].id_node = id_node;
		psorted_nodes[i].id_ver = GetCorrespondingVertex(id_node);
		psorted_nodes[i].id_cseg = GetCurrentCSeg();

		if (bClockwise)
			i++;
		else
			i = (i+num_node-1)%num_node;

		id_node = GetNextNode();
	}

	// if any hole existing
	int	id_region_hole;
	if (NumHoles()>0)
	{
		if (ppcontour_nodes_holes != NULL)
			delete []ppcontour_nodes_holes;
		ppcontour_nodes_holes = new MCMC_PSortedNodes[NumHoles()];
		assert (ppcontour_nodes_holes!=NULL);
		i = 0;

		for (int j=0; j<NumHoles(); j++)
		{
			id_region_hole = GetHole(j);

			ppcontour_nodes_holes[i++] = LayerRegion()->GetSortedNodes(id_region_hole);
		}
	}
	ComputeArea();
}

MCMC_NodeVertex* MCMC_SortedNodes::Get(const int i)
{
	int j;

	if (i>=0 && i<num_node)
		return &(psorted_nodes[i]);
	else if (i>=num_node)
	{
		int index = num_node;
		for (j=0; j<NumHoles(); j++)
		{
			int	temp = ppcontour_nodes_holes[j]->NumNodes();
			int index_sub = i-index;
			if (index_sub>=0 && index_sub<temp)
				return ppcontour_nodes_holes[j]->Get(index_sub);
			else
				index += temp;
		}
	}
	return NULL;
}

int MCMC_SortedNodes::GetFirstIndex(void)
{
	index_move = 0;
	return index_move;
}

int MCMC_SortedNodes::GetNextIndex(void)
{
	index_move++;

	MCMC_NodeVertex* ptemp = Get(index_move);
	if (ptemp!=NULL)
		return index_move;
	else
		return -1;
}

int	MCMC_SortedNodes::GetSortedIndex(const int id_node1)
{
	int i,index,index_sub;

	for (i=0; i<num_node; i++)
		if (psorted_nodes[i].id_node == id_node1)
		{
			index = i;
			return index;
		}
	// maybe in its hole
	if (i>=num_node)
	{
		index = num_node;
		for (i=0; i<NumHoles(); i++)
		{
			int	temp = ppcontour_nodes_holes[i]->NumNodes();
			index_sub = ppcontour_nodes_holes[i]->GetSortedIndex(id_node1);
			if (index_sub>=0 && index_sub<temp)
			{
				index += index_sub;
				return index;
			}
			else
				index += temp;
		}
	}
	return -1;
}

McPoint2D MCMC_SortedNodes::GetOff(const int index, McPoint2D &ofSt, McPoint2D &ofSm)
{
	assert(index>=0 && index<num_node);
	// assign values
	ofSt = psorted_nodes[index].ofStat;
	ofSm = psorted_nodes[index].ofSmooth;
	return psorted_nodes[index].ofOverall;
}

void MCMC_SortedNodes::SetOff(const int index, const McPoint2D &ofOv, const McPoint2D &ofSt,
							  const McPoint2D &ofSm)
{
	assert(index>=0 && index<num_node);
	// assign values
	psorted_nodes[index].ofStat = ofSt;
	psorted_nodes[index].ofSmooth = ofSm;
	psorted_nodes[index].ofOverall = ofOv;
}

int	MCMC_SortedNodes::GetSortedNode(const int index)
{
	assert(index>=0);
	if (index<num_node)
		return psorted_nodes[index].id_node;
	else
	{
		// on its holes
		int sub_index = index-num_node;
		for (int i=0; i<NumHoles(); i++)
		{
			int temp = ppcontour_nodes_holes[i]->NumNodes();
			if (sub_index<temp)
				return ppcontour_nodes_holes[i]->GetSortedNode(sub_index);
			else
				sub_index = sub_index-temp;
		}
	}
	return -1;
}

int	MCMC_SortedNodes::GetSortedVertex(const int index)
{
	assert(index>=0);
	if (index<num_node)
		return psorted_nodes[index].id_ver;
	else
	{
		// on its holes
		int sub_index = index-num_node;
		for (int i=0; i<NumHoles(); i++)
		{
			int temp = ppcontour_nodes_holes[i]->NumNodes();
			if (sub_index<temp)
				return ppcontour_nodes_holes[i]->GetSortedVertex(sub_index);
			else
				sub_index = sub_index-temp;
		}
	}
	return -1;
}

int	MCMC_SortedNodes::GetSortedSeg(const int index)
{
	assert(index>=0);
	if (index<num_node)
		return psorted_nodes[index].id_cseg;
	else
	{
		// on its holes
		int sub_index = index-num_node;
		for (int i=0; i<NumHoles(); i++)
		{
			int temp = ppcontour_nodes_holes[i]->NumNodes();
			if (sub_index<temp)
				return ppcontour_nodes_holes[i]->GetSortedSeg(sub_index);
			else
				sub_index = sub_index-temp;
		}
	}
	return -1;
}

int	MCMC_SortedNodes::GetSortedNextNode(const int index)
{
	assert(index>=0);
	if (index<num_node)
		return psorted_nodes[(index+1)%num_node].id_node;
	else
	{
		// on its holes
		int sub_index = index-num_node;
		for (int i=0; i<NumHoles(); i++)
		{
			int temp = ppcontour_nodes_holes[i]->NumNodes();
			if (sub_index<temp)
				return ppcontour_nodes_holes[i]->GetSortedLastNode(sub_index);
			else
				sub_index = sub_index-temp;
		}
	}
	return -1;
}

int	MCMC_SortedNodes::GetSortedLastNode(const int index)
{
	assert(index>=0);
	if (index<num_node)
		return psorted_nodes[(index+num_node-1)%num_node].id_node;
	else
	{
		// on its holes
		int sub_index = index-num_node;
		for (int i=0; i<NumHoles(); i++)
		{
			int temp = ppcontour_nodes_holes[i]->NumNodes();
			if (sub_index<temp)
				return ppcontour_nodes_holes[i]->GetSortedNextNode(sub_index);
			else
				sub_index = sub_index-temp;
		}
	}
	return -1;
}

void MCMC_SortedNodes::GiveUpData(void)
{
	num_node = 0;
	psorted_nodes = NULL;
	ppcontour_nodes_holes = NULL;
}


//////////////////////////////////////////////////////////////
// MCMC element
/////////////////////////////////////////////////////////////
MCMC_RegionElement::MCMC_RegionElement()
{
	indexid_parent	= -1;
	start_cseg		= -1;

	paffinity_map	= NULL;
	plow_level_stat	= NULL;
}

MCMC_RegionElement::~MCMC_RegionElement()
{
	Empty();
}

void MCMC_RegionElement::Empty(void)
{
	DeleteModel();
	MCMC_SortedNodes::Empty();

	dArea			= 0.0;
	indexid_parent	= -1;
	start_cseg		= -1;

	if (paffinity_map != NULL)
		delete paffinity_map;
	paffinity_map = NULL;

	if (plow_level_stat != NULL)
		delete plow_level_stat;
	plow_level_stat = NULL;
}

void MCMC_RegionElement::Init(void)
{
	if(paffinity_map == NULL)
		paffinity_map = new MCMC_AffinityMap;

	if(plow_level_stat == NULL)
		plow_level_stat = new MCMC_LowLevelStat4Region;
}


MCMC_RegionElement& MCMC_RegionElement::operator=(const MCMC_RegionElement& a)
{
	if (this != &a)
	{
		*((MCMC_Region_Stat*) this) = *((MCMC_Region_Stat*) &a);

		start_cseg		= a.start_cseg;
		indexid_parent	= a.indexid_parent;

		direction		= a.direction;
		clockwise		= a.clockwise;
		last_seg		= a.last_seg;
		current_ver		= a.current_ver;
		id_node_temp	= a.id_node_temp;
		id_node_first	= a.id_node_first;
		energy_shapeprior = a.energy_shapeprior;
		pt_center		= a.pt_center;

		if (paffinity_map != NULL)
			delete paffinity_map;
		paffinity_map = NULL;
		if (a.paffinity_map != NULL)
		{
			paffinity_map = new MCMC_AffinityMap;
			*paffinity_map = *a.paffinity_map;
		}

		if (plow_level_stat != NULL)
			delete plow_level_stat;
		plow_level_stat = NULL;
		if (a.plow_level_stat != NULL)
		{
			plow_level_stat = new MCMC_LowLevelStat4Region;
			*plow_level_stat = *a.plow_level_stat;
		}
	}

	return *this;
}

ofstream & operator>>(MCMC_RegionElement &a, ofstream &strm)
{
	*((MCMC_Region_Stat*)&a) >> strm;

	strm.write((char*)&(a.start_cseg), sizeof(a.start_cseg));
	strm.write((char*)&(a.indexid_parent), sizeof(a.indexid_parent));

	return strm;
}

ifstream & operator<<(MCMC_RegionElement &a, ifstream &strm)
{
	*((MCMC_Region_Stat*)&a) << strm;

	strm.read((char*)&(a.start_cseg), sizeof(a.start_cseg));
	strm.read((char*)&(a.indexid_parent), sizeof(a.indexid_parent));

	return strm;
}

// get current contour segment
int	MCMC_RegionElement::GetCurrentCSeg(void)
{
	if (!Valid())
		return -1;

	if (current_ver >= 0)
	{
		int current_seg;
		
		// get current vertex pointer
		MCMC_Vertex* pvertex = LayerRegion()->GetVertex(current_ver);
		mcmcVector<int> vc_csegments;

		GetCSetFromVertex(pvertex, vc_csegments, IndexId());
		if (last_seg==-1)
			current_seg = start_cseg;
		else
			current_seg = WhichWay(last_seg, vc_csegments);
		return current_seg;
	}
	else
		return -1;
}


int	MCMC_RegionElement::GetFirstCSeg(void)
{
	MCMC_Contour_Segment *pcseg;

	pcseg = LayerRegion()->GetCSeg(start_cseg);
	if (pcseg == NULL)
		return -1;

	if (!IsTheSameElement(pcseg->GetLeftAtomicRegion(), IndexId()) && 
		!IsTheSameElement(pcseg->GetRightAtomicRegion(), IndexId()))
		FindTheStartCSeg();
	
	if (GetFirstVer()<0)
		return -1;
	else
		return GetCurrentCSeg();
}
// get the next contour segment of this region
int	MCMC_RegionElement::GetNextCSeg(void)
{
	if (GetNextVer()<0)
		return -1;
	else
		return GetCurrentCSeg();
}

// get the first vertex of this region
int	MCMC_RegionElement::GetFirstVer(void)
{
	if (!Valid())
		return -1;

	int	current_seg = start_cseg;

	MCMC_Contour_Segment* pcseg = LayerRegion()->GetCSeg(current_seg);
	assert(pcseg != NULL);
	if (pcseg==NULL)
		return -1;
	int iVer1,iVer2;
	pcseg->GetVertices(iVer1,iVer2);
	if (IsTheSameElement(pcseg->GetRightAtomicRegion(), IndexId()))
		current_ver = iVer1;
	else
		current_ver = iVer2;

	last_seg = -1;
	return	current_ver;
}


// get the next vertex of this region
int	MCMC_RegionElement::GetNextVer(void)
{
	if (!Valid())
		return -1;

	int	current_seg;

	if (current_ver>=0)
	{
		// get current vertex
		MCMC_Vertex*	pvertex = LayerRegion()->GetVertex(current_ver);
		mcmcVector<int>		vc_csegments;

		GetCSetFromVertex(pvertex, vc_csegments, IndexId());

		if (last_seg==-1)
			current_seg = start_cseg;
		else
			current_seg = WhichWay(last_seg, vc_csegments);
		if (current_seg<0)
			return -1;
		last_seg = current_seg;
		// get corresponding contour segment
		MCMC_Contour_Segment* pcseg = LayerRegion()->GetCSeg(current_seg);
		assert(pcseg != NULL);
		// go to the next vertex
		current_ver = pcseg->GetTheOtherVertex(current_ver);
		pvertex = LayerRegion()->GetVertex(current_ver);

		// get corresponding segment
		GetCSetFromVertex(pvertex, vc_csegments, IndexId());
		current_seg = WhichWay(last_seg, vc_csegments);
		if (current_seg==start_cseg || current_seg<0)	// go back to starting point
			return -1;
		else
			return current_ver;
	}
	return -1;
}


int MCMC_RegionElement::ToppestLevelIndexId(void)
{
	MCMC_RegionElement *pregion;

	pregion = this;
	while (pregion->IndexIdParent()>=0)
	{
		pregion = LayerRegion()->GetRegionFromIndexId(pregion->IndexIdParent());
		assert(pregion != NULL);
	}

	return pregion->IndexId();
}


bool MCMC_RegionElement::IsInGRegion(const int indexid_gregion)
{
	// valid id_gregion?
	if (indexid_gregion<0 || LayerRegion()==NULL || indexid_parent<0)
		return false;
	// id
	if (IndexId()==indexid_gregion || indexid_parent==indexid_gregion)
		return true;

	MCMC_RegionElement	*pregion_element;

	pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_parent);
	while (pregion_element != NULL)
	{
		if (pregion_element->IndexIdParent() == indexid_gregion)
			return true;

		pregion_element = LayerRegion()->GetRegionFromIndexId(pregion_element->IndexIdParent());
	}
	return false;
}

int	MCMC_RegionElement::WhichWay(int cseg_last, mcmcVector<int> &vc_csegments)
{
	if (vc_csegments.size()==1 || cseg_last==-1)
		return vc_csegments(0);		// no other choices
	else
	{
		MCMC_Contour_Segment *pcseg1,*pcseg2;
		double	x1,y1,x2,y2,x3,y3;
		int		id_node1,id_node2,id_node3;

		pcseg1 = LayerRegion()->GetCSeg(cseg_last);
		assert(pcseg1 != NULL);

		if (IsTheSameElement(pcseg1->GetLeftAtomicRegion(), IndexId()))
			pcseg1->SwapVertices();
		id_node1 = pcseg1->extn[1];
		id_node2 = pcseg1->node[1];

		MCMC_Node	*pnode;
		double		angle1,angle2,angle_temp=1e10;
		int			itemp=-1;

		pnode = LayerRegion()->GetNode(id_node1);
		assert(pnode != NULL);
		x1 = pnode->GetX();
		y1 = pnode->GetY();
		pnode = LayerRegion()->GetNode(id_node2);
		assert(pnode != NULL);
		x2 = pnode->GetX();
		y2 = pnode->GetY();
		angle1 = xysa(x2,y2,x1,y1);

		for (int i=0; i<vc_csegments.size(); i++)
		{
			pcseg2 = LayerRegion()->GetCSeg(vc_csegments(i));
			if(!(pcseg2->SelfClosed() && pcseg2->Id()!=cseg_last) &&
			   pcseg2->GetLeftAtomicRegion()!=pcseg2->GetRightAtomicRegion())
			{
				if (IsTheSameElement(pcseg2->GetLeftAtomicRegion(), IndexId()))
					pcseg2->SwapVertices();
				id_node3 = pcseg2->extn[0];
				pnode = LayerRegion()->GetNode(id_node3);
				assert(pnode != NULL);
				x3 = pnode->GetX();
				y3 = pnode->GetY();
				angle2 = xysa(x2,y2,x3,y3);

				assert(pcseg2 != NULL);
				if (AngleO2PI(angle1-angle2)<angle_temp	)
				{
					itemp = i;
					angle_temp = AngleO2PI(angle1-angle2);
				}
			}
		}
		if (itemp>=0)
			return vc_csegments(itemp);
		else
			return -1;
	}
}

int	MCMC_RegionElement::GetFirstNode()
{
	int id_cseg;

	id_cseg = GetFirstCSeg();
	if (id_cseg < 0)
		return -1;
	else
	{
		MCMC_Contour_Segment	*pcseg;
		int						id_node;

		// get current contour segment
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);

		// get the first node on the segment
		if (pcseg->end[0] == current_ver)
		{
			id_node = pcseg->GetFirstNode();
			direction = 1;
		}
		else
		{
			id_node = pcseg->GetLastNode();
			direction = 0;
		}

		id_node_temp = id_node;
		id_node_first = id_node;
		return id_node;
	}
}

int	MCMC_RegionElement::GetNextNode(void)
{
	// get current segment
	int id_cseg = GetCurrentCSeg();
	if (id_cseg < 0)
		return -1;
	else
	{
		int					 id_node;
		MCMC_Contour_Segment *pcseg;

		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);

		if (direction)
			id_node = pcseg->GetNextNode();
		else
			id_node = pcseg->GetPrevNode();
		if (id_node >= 0)
		{
			id_node_temp = id_node;
			if (id_node != id_node_first)
				return id_node;		// get the next node
			else
				return -1;
		}
		else
		{
			// start from the next segment
			id_cseg = GetNextCSeg();
			if (id_cseg < 0)
				return -1;
			else
			{
				pcseg = LayerRegion()->GetCSeg(id_cseg);
				assert(pcseg != NULL);

				// skip the starting node
				id_node = pcseg->GetFirstNode();
				if (id_node == id_node_temp)
					direction = 1;
				else
				{
					id_node = pcseg->GetLastNode();
					assert (id_node == id_node_temp);
					direction = 0;
				}
				if (direction)
					id_node = pcseg->GetNextNode();
				else
					id_node = pcseg->GetPrevNode();

				id_node_temp = id_node;
				if (id_node != id_node_first)
					return id_node;
				else
					return -1;
			}
		}
	}
}

int	MCMC_RegionElement::NumNodes(void)
{
	int id_node,i=0,max_node;
	
	max_node = LayerRegion()->NumNodes();
	
	id_node = GetFirstNode();
	while (id_node >= 0)
	{
		i++;
		
		if (i > max_node)
		{
			cout << "NumNodes" << endl;
			XPT_OutOfBound  *pxpt;

			pxpt = new XPT_OutOfBound("topology error");
			throw pxpt;
		}

		id_node = GetNextNode();
	}
	return i;
}

static int grid_index_left_row[3][3] ={{0,-1,0},{0, 0,-1},{0, 0,0}};
static int grid_index_left_col[3][3] ={{0,-1,0},{-1,0, 0},{0, 0,0}};
static int grid_index_right_row[3][3]={{0,-1,0},{-1,0, 0},{0, 0,0}};
static int grid_index_right_col[3][3]={{0, 0,0},{-1,0, 0},{0,-1,0}};

inline void	MCMC_RegionElement::GridLeftRightRegionId(int &id_region_left, McImagePoint &pt_left_pixel,
									  int &id_region_right, McImagePoint &pt_right_pixel,
									  const int row1, const int col1, const int row2, const int col2)
{
	int row_region,col_region;
	int index_row,index_col;

	index_row = row2-row1+1;
	index_col = col2-col1+1;

	// left region
	row_region = row1+grid_index_left_row[index_row][index_col];
	col_region = col1+grid_index_left_col[index_row][index_col];
	if (player_region->mx_label_id_aregions.Valid(row_region, col_region))
		id_region_left = player_region->mx_label_id_aregions(row_region, col_region);
	else
		id_region_left = -1;
	pt_left_pixel.SetRow(row_region);
	pt_left_pixel.SetCol(col_region);

	// right region
	row_region = row1+grid_index_right_row[index_row][index_col];
	col_region = col1+grid_index_right_col[index_row][index_col];
	if (player_region->mx_label_id_aregions.Valid(row_region, col_region))
		id_region_right = LayerRegion()->mx_label_id_aregions(row_region, col_region);
	else
		id_region_right = -1;

	pt_right_pixel.SetRow(row_region);
	pt_right_pixel.SetCol(col_region);
}

MCMC_Node* MCMC_RegionElement::GridGetNextNode(MCMC_Node *pnode1, McImagePoint &pt_left_pixel, McImagePoint &pt_right_pixel, const bool bclockwise)
{
	assert(pnode1 != NULL);

	MCMC_Node	*pnode2;
	static int	p_id_nodes[2];
	int			i,id_left_region,id_right_region,id_node2;


	// check the link first
	p_id_nodes[0] = pnode1->link[0];
	p_id_nodes[1] = pnode1->link[1];
	for (i=0; i<2; i++)
	{
		id_node2 = p_id_nodes[i];
		if (id_node2 >= 0)
		{
			pnode2 = player_region->GetNode(id_node2);
			assert(pnode2 != NULL);
			if (abs(pnode1->GridRow()-pnode2->GridRow())+abs(pnode1->GridCol()-pnode2->GridCol()) == 1)
			{
				// the neighboring nodes on the grid
				GridLeftRightRegionId(id_left_region, pt_left_pixel, id_right_region, pt_right_pixel,
									  pnode1->GridRow(), pnode1->GridCol(), pnode2->GridRow(), pnode2->GridCol());
				if (bclockwise && IsTheSameElement(id_right_region, IndexId()) && !IsTheSameElement(id_left_region, IndexId()))
					return pnode2;
				if (!bclockwise && IsTheSameElement(id_left_region, IndexId()) && !IsTheSameElement(id_right_region, IndexId()))
					return pnode2;
			}
		}
	}

	// search the neighbors
	int row1,col1,row2,col2;

	row1 = pnode1->GridRow();
	col1 = pnode1->GridCol();

	for (i=0; i<4; i++)
	{
		row2 = row1+row_neighbors4[i];
		col2 = col1+col_neighbors4[i];

		if (player_region->mx_label_id_nodes.Valid(row2, col2))
		{
			pnode2 = player_region->GetNode(player_region->mx_label_id_nodes(row2, col2));

			// the neighboring nodes on the grid
			GridLeftRightRegionId(id_left_region, pt_left_pixel, id_right_region, pt_right_pixel, row1, col1, row2, col2);
			if (bclockwise && IsTheSameElement(id_right_region, IndexId()) && !IsTheSameElement(id_left_region, IndexId()))
				return pnode2;
			if (!bclockwise && IsTheSameElement(id_left_region, IndexId()) && !IsTheSameElement(id_right_region, IndexId()))
				return pnode2;
		}
	}

	return NULL;
}

void MCMC_RegionElement::GridGetNeighborNodes(MCMC_Node *pnode1, MCMC_Node *p_nodes_left[], MCMC_Node *p_nodes_right[], const int num)
{
	assert(pnode1 != NULL);

	MCMC_Node		*pnode_temp;
	int				i;
	McImagePoint	pt_left_pixel,pt_right_pixel;

	pnode_temp = pnode1;
	for (i=0; i<num; i++)
	{
		// clockwise on the right side
		p_nodes_right[i] = GridGetNextNode(pnode_temp, pt_left_pixel, pt_right_pixel, true);
		//		assert_ztu(p_nodes_right[i] != NULL);

		if (p_nodes_right[i] == NULL)
		{
			p_nodes_right[0] = NULL;
			return;
		}


		pnode_temp = p_nodes_right[i];
	}

	pnode_temp = pnode1;
	for (i=0; i<num; i++)
	{
		// clockwise on the right side
		p_nodes_left[i] = GridGetNextNode(pnode_temp, pt_left_pixel, pt_right_pixel, false);
		//		assert_ztu(p_nodes_left[i] != NULL);

		if (p_nodes_left[i] == NULL)
		{
			p_nodes_left[0] = NULL;
			return;
		}


		pnode_temp = p_nodes_left[i];
	}
}

double MCMC_RegionElement::GridGetStatForce(MCMC_Node *pnode, const McImagePoint pt_pixel)
{
	if (Model() != NULL)
		return Model()->GetLogPostProb(pnode->GridRow(), pnode->GridCol());
	else
	{
		MCMC_RegionElement *pregion;
		
		pregion = LayerRegion()->GetRegionFromIndexId(IndexIdParent());
		while (pregion != NULL)
		{
			if (pregion->Model() != NULL)
				return pregion->Model()->GetLogPostProb(pnode->GridRow(), pnode->GridCol());

			pregion = LayerRegion()->GetRegionFromIndexId(pregion->IndexIdParent());
		}
		return 0.0;
	}
}

double MCMC_RegionElement::GridGetShapeForce(MCMC_Node *pnode, MCMC_Node *p_nodes_left[], MCMC_Node *p_nodes_right[])
{
	double		curvature, norm_direction,grid_direction;
	McPoint2D	pt1,pt,pt2;

	// simple smoothness force
	// the three points
	pt1.Set(p_nodes_left[1]->GetX(), p_nodes_left[1]->GetY());
	pt.Set(pnode->GetX(), pnode->GetY());
	pt2.Set(p_nodes_right[1]->GetX(), p_nodes_right[1]->GetY());

	norm_direction = ComputeNormalCurvature(pt1, pt, pt2, &curvature);
	grid_direction = atan2(p_nodes_right[0]->GetY()-pnode->GetY(), p_nodes_right[0]->GetX()-pnode->GetX())+PI/2.0;

	return curvature*(cos(norm_direction)*cos(grid_direction)+sin(norm_direction)*sin(grid_direction));
}

double MCMC_RegionElement::GridGetTotalForce(MCMC_Node *pnode, const McImagePoint pt_pixel, MCMC_Node *p_nodes_left[], MCMC_Node *p_nodes_right[])
{
	double dforce_stat,dforce_shape,dforce_total;

	dforce_stat = GridGetStatForce(pnode, pt_pixel);
	dforce_shape = GridGetShapeForce(pnode, p_nodes_left, p_nodes_right);
	
	dforce_total = (dforce_stat+0.5*MU*dforce_shape)/FORCE_NORMAL;

	return dforce_total;
}

double MCMC_RegionElement::GridGetDiffusionMove(const int id_node1, int &id_node2)
{
	MCMC_Node		*pnode1,*pnode2,*p_nodes_right[3],*p_nodes_left[3];
	McImagePoint	pt_left_pixel,pt_right_pixel;
	double			dforce_this_region,dforce_neighbor_region;
	MCMC_AtomicRegion *paregion_neighbor;

	id_node2 = -1;

	pnode1 = player_region->GetNode(id_node1);
	if (pnode1 == NULL)
		return 0.0;
	// the next node on the grid
	pnode2 = GridGetNextNode(pnode1, pt_left_pixel, pt_right_pixel, true);
	if (pnode2 != NULL)
		id_node2 = player_region->mx_label_id_nodes(pnode2->GetRow(), pnode2->GetCol());

	if (pnode2 == NULL || !player_region->mx_label_id_aregions.Valid(pt_right_pixel.Row(), pt_right_pixel.Col()) ||
						  !player_region->mx_label_id_aregions.Valid(pt_left_pixel.Row(), pt_left_pixel.Col()))
		return 0.0;

	paregion_neighbor = player_region->GetAtomicRegion(player_region->mx_label_id_aregions(pt_left_pixel.Row(), pt_left_pixel.Col()));
	if (paregion_neighbor == NULL)
		return 0.0;

	// neighbor nodes
	GridGetNeighborNodes(pnode1, p_nodes_left, p_nodes_right, 3);
	if (p_nodes_left[0] == NULL || p_nodes_right[0] == NULL)
	{
		id_node2 = -1;
		XPT_OutOfBound  *pxpt;
		pxpt = new XPT_OutOfBound("diffusion problem");

		throw pxpt;
	}

	dforce_this_region = GridGetTotalForce(pnode1, pt_left_pixel, p_nodes_left, p_nodes_right);

	// neighbor nodes
	paregion_neighbor->GridGetNeighborNodes(pnode2, p_nodes_left, p_nodes_right, 3);
	if (p_nodes_left[0] == NULL || p_nodes_right[0] == NULL)
	{
		id_node2 = -1;
		XPT_OutOfBound  *pxpt;
		pxpt = new XPT_OutOfBound("diffusion problem");

		throw pxpt;
	}

	dforce_neighbor_region = paregion_neighbor->GridGetTotalForce(pnode2, pt_right_pixel, p_nodes_left, p_nodes_right);

	return dforce_this_region-dforce_neighbor_region;
}


#define MAX_DIST_TO_CONSIDER  30

double	MCMC_RegionElement::area(SimpleMatrix<McPoint2D> &vc_points)
{
	double	dResult = 0;
	int		i,cnt;
	
	cnt = vc_points.size();
	if (cnt > 2)
	{
		// valid region
		for (i = 0; i < cnt; i++)
			dResult += vc_points(i).x*(vc_points((i+cnt-1)%cnt).y-vc_points((i+1)%cnt).y);
		dResult = fabs(dResult)/2;
	}
	return dResult;
}

double MCMC_RegionElement::MinDistToPoints(SimpleMatrix<McPoint2D> &vc_points, const int i)
{
	// to compute the normal direction
	double	angle1,angle2,angle_normal;
	int		j;
	
	// the previous point;
	j = (i-2+vc_points.size())%vc_points.size();
	angle1 = atan2(vc_points(j).y-vc_points(i).y, vc_points(j).x-vc_points(i).x);

	j = (i+2)%vc_points.size();
	angle2 = atan2(vc_points(j).y-vc_points(i).y, vc_points(j).x-vc_points(i).x);

	angle_normal = angle1+AngleO2PI(angle2-angle1)/2.0;
	
	McPoint2D	pt11,pt12,pt21,pt22,pt_cross;
	double		dist_min,dist_temp;

	dist_min = MAX_DIST_TO_CONSIDER;

	// current point
	pt11 = vc_points(i);
	// to cross
	pt12.x = pt11.x+MAX_DIST_TO_CONSIDER*cos(angle_normal);
	pt12.y = pt11.y+MAX_DIST_TO_CONSIDER*sin(angle_normal);

	for (j=(i+2)%vc_points.size(); j!=(i-2+vc_points.size())%vc_points.size(); j=(j+1)%vc_points.size())
	{
		pt21 = vc_points(j);
		if (Distance(pt11, pt21) < MAX_DIST_TO_CONSIDER)
		{
			pt22 = vc_points((j+1)%vc_points.size());

			if (LineIntersection(&pt_cross, pt11, pt12, pt21, pt22) == CROSSING)
			{
				dist_temp = Distance(pt_cross, pt11);
				if (dist_temp < dist_min)
					dist_min = dist_temp;
			}

		}
	}
	
	return dist_min;
}

double MCMC_RegionElement::ComputeEnergyShape(void)
{
	energy_shapeprior = GetLength()*PRIOR_SHAPE + PRIOR_REGION;

	return energy_shapeprior;
}


void MCMC_RegionElement::GetRegionDomain(LinkList<PMCMC_PixelSet> &list_domain)
{
	// initialization
	list_domain.Empty();

	MCMC_AtomicRegion *paregion;

	paregion = LayerRegion()->GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		list_domain.Insert(&paregion->pixel_set);

		paregion = LayerRegion()->GetAtomicRegion(GetNextAtomicRegion());
	}

}

void MCMC_RegionElement::LabelFromDomain(Raster<int> &ras_label_map, int id_label)
{
	bool	bgo;
	int		i,j;

	bgo = GetFirstPixel(j, i);
	while (bgo)
	{
		ras_label_map(j, i) = id_label;
		bgo = GetNextPixel(j, i);
	}
}


//////////////////////////////////////////////////////////////
// MCMC Region
/////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////
// construction & destruction
MCMC_AtomicRegion::MCMC_AtomicRegion()
{
	energy_shapeprior = 0.0;
}

MCMC_AtomicRegion::~MCMC_AtomicRegion()
{
	Empty();
}


void MCMC_AtomicRegion::Init(void)
{
	// free allocated space
	Empty();

	MCMC_RegionElement::Init();

	// holes
	holes.Init();

	// neighbors
	neighbors.Init();
}


void MCMC_AtomicRegion::Empty(void)
{
	MCMC_RegionElement::Empty();

	current_ver = -1;		// current segment of this region
	start_cseg = -1;

	// holes
	holes.Free();

	// neighbors
	neighbors.Free();

	energy_shapeprior = 0.0;
}

MCMC_AtomicRegion& MCMC_AtomicRegion::operator=(MCMC_AtomicRegion& a)
{
	if (this != &a)
	{
		*((MCMC_RegionElement*) this) = *((MCMC_RegionElement*) &a);

	
		pixel_set		= a.pixel_set;
		holes			= a.holes;
		neighbors		= a.neighbors;
		alive			= a.alive;
	}

	return *this;
}

int MCMC_AtomicRegion::TotalCodingBits(void)
{
	int nbytes;

	nbytes = sizeof(start_cseg)+
			sizeof(clockwise)+
			sizeof(last_seg);

	if (Model()!=NULL)
		nbytes += Model()->TotalCodingBits();
	return nbytes;
}


ofstream & operator>>(MCMC_AtomicRegion &a, ofstream &strm)
{
	*((MCMC_RegionElement*)&a) >> strm;
	*((MCMC_RegionElement*)&a)>> strm;

	// holes
	a.holes >> strm;
	// neighbors
	a.neighbors >> strm;

	strm.write((char*)&(a.alive), sizeof(a.alive));
	strm.write((char*)&(a.energy_shapeprior), sizeof(a.energy_shapeprior));

	return strm;
}

ifstream & operator<<(MCMC_AtomicRegion &a, ifstream &strm)
{
	a.Init();

	*((MCMC_RegionElement*)&a) << strm;
	*((MCMC_RegionElement*)&a) << strm;

	// holes
	a.holes << strm;
	// neighbors
	a.neighbors << strm;

	strm.read((char*)&(a.alive), sizeof(a.alive));
	strm.read((char*)&(a.energy_shapeprior), sizeof(a.energy_shapeprior));

	return strm;
}


//////////////////////////////////////////////////////////////
// access

void MCMC_AtomicRegion::FindTheStartCSeg(void)
{
	int id_cseg_buf;
	MCMC_Contour_Segment *pcseg;

	id_cseg_buf = LayerRegion()->management_csegments.current;
	pcseg = LayerRegion()->GetCSeg(LayerRegion()->GetFirstCSeg());
	while (pcseg != NULL)
	{
		if (pcseg->BelongToAtomicRegion(Id()))
		{
			ChangeStartCSeg(pcseg->Id());
			break;
		}
		pcseg = LayerRegion()->GetCSeg(LayerRegion()->GetNextCSeg());
	}
	if (pcseg == NULL)				
		ChangeStartCSeg(-1);

	LayerRegion()->management_csegments.current = id_cseg_buf;
}


// side 0: right  1:left)
void MCMC_AtomicRegion::GetCSetFromVertex(MCMC_Vertex *pver, mcmcVector<int> &vc_csegs, const int id_region, const int side)
{
	if (pver != NULL)
		pver->GetCSegOfAtomicRegion(vc_csegs, id_region, side);
}


// Is one of its holes?
bool MCMC_AtomicRegion::IsHole(const int id_rgn)
{
	if (id_rgn<0 || id_rgn==Id())
		return false;

	return holes.Exist(id_rgn);
}

// add one region to this region as a hole
bool MCMC_AtomicRegion::AddHole(const int id_rgn)
{
	if (id_rgn<0 || id_rgn==Id())
		return false;

	return holes.Add(id_rgn);
}

// Is one its neighbors?
bool MCMC_AtomicRegion::IsNeighbor(const int id_rgn)
{
	if (id_rgn<0 || id_rgn==Id())
		return false;

	return neighbors.Exist(id_rgn);
}

// add one region to this region as a neighbor
bool MCMC_AtomicRegion::AddNeighbor(const int id_rgn)
{
	if (id_rgn<0 || id_rgn==Id())
		return false;

	return neighbors.Add(id_rgn);
}

bool MCMC_AtomicRegion::DelNeighbor(const int id_region_nb,const int id_region_replace)
{
	MCMC_AtomicRegion *pregion_replace;

	pregion_replace = LayerRegion()->GetAtomicRegion(id_region_replace);

	if (pregion_replace != NULL)
		pregion_replace->AddNeighbor(Id());

	if (id_region_replace == Id())
		return neighbors.Delete(id_region_nb);
	else
		return neighbors.Delete(id_region_nb, id_region_replace);
}

void MCMC_AtomicRegion::DelFromNb(const int id_region_replace)
{
	int			id_region_nb;
	MCMC_AtomicRegion	*pregion_nb;

	id_region_nb = GetFirstNeighbor();
	while (id_region_nb>=0)
	{
		pregion_nb = LayerRegion()->GetAtomicRegion(id_region_nb);
		assert(pregion_nb != NULL);
		pregion_nb->DelNeighbor(Id(),id_region_replace);

		id_region_nb = GetNextNeighbor();
	}
}

// rearrange segments, merge segments that have one degree vertex
void MCMC_AtomicRegion::MergeSegmentsOfRegion(void)
{
	int			id_cseg_start,id_cseg1,id_cseg2;

	if (!Valid())
		return;

	id_cseg_start = GetFirstCSeg();
	id_cseg1 = id_cseg_start;
	id_cseg2 = GetNextCSeg();
	while (id_cseg1>=0 && id_cseg2>=0)
	{
		if (LayerRegion()->MergeSegments(id_cseg1,id_cseg2,Id())==1)
		{
			id_cseg_start = GetFirstCSeg();
			id_cseg1 = id_cseg_start;
		}
		else
			id_cseg1 = id_cseg2;
		id_cseg2 = GetNextCSeg();
	}
	LayerRegion()->MergeSegments(id_cseg1,id_cseg_start,Id());
}


int MCMC_AtomicRegion::Label(Raster<int> &ras_label_map, int id_label)
{
	MCMC_Block block;

	return block.Label(*this,ras_label_map,id_label);
}

// create the pxiel set for this region
void MCMC_AtomicRegion::CreateRegionDomain(void)
{
	MCMC_Block block;

	pixel_set.Empty();
	block.LabelToPixelSet(*this, pixel_set);

	bool	bgo;
	int		i,j;

	pt_center.Set(0, 0);
	bgo = pixel_set.GetFirst(j, i);
	while (bgo)
	{
		pt_center.x += j;
		pt_center.y += i;

		bgo = pixel_set.GetNext(j, i);
	}
	if (pixel_set.Num() > 0)
		pt_center = pt_center/((double)pixel_set.Num());
}


int	MCMC_AtomicRegion::GetVerNum(void)
{
	int i=0,id_ver;

	id_ver = GetFirstVer();
	while (id_ver >= 0)
	{
		i++;
		id_ver = GetNextVer();
	}
	return i;
}

bool MCMC_AtomicRegion::SegOfHole(const int id_cseg)
{
	bool	bResult=false;
	int		i;
	MCMC_Contour_Segment *pcseg;


	pcseg = LayerRegion()->GetCSeg(id_cseg);
	if (pcseg != NULL)
	{
		int id_nb;


		if (pcseg->boundry_of[0] == Id())
			id_nb = pcseg->boundry_of[1];
		else
		{
			assert(pcseg->boundry_of[1] == Id());
			id_nb = pcseg->boundry_of[0];
		}

		for (i=0; i<holes.Num(); i++)
		{
			if (holes.Data(i) == id_nb)
			{
				bResult = true;
				break;
			}
		}
	}
	return bResult;
}


// flag=0, set to -1 if no other contour segment/*
void MCMC_AtomicRegion::UseDifferentStartCSeg(void)
{
	GetFirstCSeg();

	MCMC_Contour_Segment *pcseg;
	pcseg = LayerRegion()->GetCSeg(GetNextCSeg());
	if (pcseg!=NULL && pcseg->Id()!=GetStartCSeg() && pcseg->BelongToAtomicRegion(Id()))
		ChangeStartCSeg(pcseg->Id());
	else
	{
		int id_cseg_buf;
		MCMC_Contour_Segment *pcseg;

		id_cseg_buf = LayerRegion()->management_csegments.current;
		pcseg = LayerRegion()->GetCSeg(LayerRegion()->GetFirstCSeg());
		while (pcseg != NULL)
		{
			if (pcseg->BelongToAtomicRegion(Id()) && pcseg->Id()!=GetStartCSeg())
			{
				ChangeStartCSeg(pcseg->Id());
				break;
			}
			pcseg = LayerRegion()->GetCSeg(LayerRegion()->GetNextCSeg());
		}
		if (pcseg == NULL)				
			ChangeStartCSeg(-1);

		LayerRegion()->management_csegments.current = id_cseg_buf;
	}
}


// get the length of its outer boundary
double MCMC_AtomicRegion::GetBoundLength(void)
{
	int					id_cseg;
	MCMC_Contour_Segment *pcseg;
	double				dLength=0.0;

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);

		dLength += pcseg->GetLength();
		id_cseg = GetNextCSeg();
	}
	return dLength;
}

// get the length of its total boundary
double MCMC_AtomicRegion::GetLength(void)
{
	double	dLength=GetBoundLength();
	int		id_region;
	MCMC_AtomicRegion *pregion;

	
	id_region = GetFirstHole();
	while (id_region >= 0)
	{
		pregion = LayerRegion()->GetAtomicRegion(id_region);
		assert(pregion != NULL);
		dLength += pregion->GetBoundLength();

		id_region = GetNextHole();
	}
	
	return dLength;
}

void MCMC_AtomicRegion::ChangeForDiffusion(void)
{
	int id_cseg;
	MCMC_Contour_Segment *pcseg;

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		if (pcseg->boundry_of[0]==-1)
			pcseg->boundry_of[0] = 0;
		if (pcseg->boundry_of[1]==-1)
			pcseg->boundry_of[1] = 0;
		id_cseg = GetNextCSeg();
	}
}

void MCMC_AtomicRegion::ChangeBack(void)
{
	int id_cseg;
	MCMC_Contour_Segment *pcseg;

	id_cseg = GetFirstCSeg();
	while (id_cseg>=0)
	{
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		if (pcseg->boundry_of[0]==0)
			pcseg->boundry_of[0] = -1;
		if (pcseg->boundry_of[1]==0)
			pcseg->boundry_of[1] = -1;
		id_cseg = GetNextCSeg();
	}
}

void MCMC_AtomicRegion::CopyRegionInfo(MCMC_AtomicRegion &region, int *pids)
{
	// copy statistics info.
	GetStat() = region.GetStat();
}

void MCMC_AtomicRegion::Smooth(void)
{
	int						id_cseg;
	MCMC_Contour_Segment	*pcseg;

	if (Valid())
	{
		id_cseg = GetFirstCSeg();
		while (id_cseg >= 0)
		{
			pcseg = LayerRegion()->GetCSeg(id_cseg);
			assert(pcseg != NULL);
			pcseg->Smooth();

			id_cseg = GetNextCSeg();
		}
	}
}

void MCMC_AtomicRegion::SetDirection(void)
{
	int						id_ver1,id_ver2,id_cseg;
	MCMC_Vertex				*pver1,*pver2;
	MCMC_Contour_Segment	*pcseg;

	id_ver1 = GetFirstVer();
	pver1 = LayerRegion()->GetVertex(id_ver1);
	if (pver1 == NULL)
		return;

	id_ver2 = GetNextVer();
	if (id_ver2<0)
		id_ver2 = id_ver1;		// one closed contour segment
	pver2 = LayerRegion()->GetVertex(id_ver2);
	assert(pver2 != NULL);

	id_cseg = GetFirstCSeg();
	pcseg = LayerRegion()->GetCSeg(id_cseg);
	assert(pcseg != NULL);

	if (! ((pcseg->end[0]==id_ver1 && pcseg->end[1]==id_ver2) ||
			(pcseg->end[0]==id_ver2 && pcseg->end[1]==id_ver1)))
	{
		cout << "SetDirection" << endl;
		XPT_OutOfBound *pxpt;
		pxpt = new XPT_OutOfBound("segment id");

		throw pxpt;
	}
	if (id_ver1 != id_ver2)
	{
		if (pcseg->end[0]==id_ver1 && pcseg->end[1]==id_ver2)
		{
			if (Id()==pcseg->GetRightAtomicRegion())
				clockwise = true;
			else
				clockwise = false;
		}
		else
		{
			if (Id()==pcseg->GetRightAtomicRegion())
				clockwise = false;
			else
				clockwise = true;
		}
	}
	else
	{
		if (Id()==pcseg->GetRightAtomicRegion())
			clockwise = true;
		else
			clockwise = false;
	}
}

double MCMC_AtomicRegion::GetNormal(const int id_node1, const int id_node, const int id_node2, double &curv, double *pmin_dis)
{
	MCMC_Node	*pnode,*pnode1,*pnode2;
	double		d1,d2,norm_direction;

	// initialization
	norm_direction = 0.0;
	curv = 0.0;
	pnode = LayerRegion()->GetNode(id_node);
	if (pnode==NULL)
		return curv;
	pnode1 = LayerRegion()->GetNode(id_node1);
	assert(pnode1 != NULL);
	pnode2 = LayerRegion()->GetNode(id_node2);
	assert(pnode2 != NULL);

	McPoint2D pt1,pt,pt2;

	pt1.x	= pnode1->GetX();		pt1.y	= pnode1->GetY();
	pt.x	= pnode->GetX();		pt.y	= pnode->GetY();
	pt2.x	= pnode2->GetX();		pt2.y	= pnode2->GetY();

	// compute the normal and curvature
	norm_direction = ComputeNormalCurvature(pt1, pt, pt2, &curv);

	// compute the curvature
	d1 = Distance(pnode->GetX(),pnode->GetY(),pnode1->GetX(),pnode1->GetY());
	d2 = Distance(pnode->GetX(),pnode->GetY(),pnode2->GetX(),pnode2->GetY());

	if (pmin_dis != NULL)
		*pmin_dis = min(d1,d2);

	return norm_direction;
}

double MCMC_AtomicRegion::GetNormal(const int index,double &curv,double *pmin_dis)
{
	int id_node1,id_node,id_node2;

	// make sure id_node_cl1, id_node and id_node2_cl appear in clockwise order
	id_node1 = GetSortedLastNode(index);
	id_node = GetSortedNode(index);
	id_node2 = GetSortedNextNode(index);

	return GetNormal(id_node1,id_node,id_node2,curv,pmin_dis);
}

double MCMC_AtomicRegion::GetStatForce(const int id_node)
{
	double			dForce;
	int				i,j,x,y,num,x1,y1;
	MCMC_Node		*pnode;

	pnode = LayerRegion()->GetNode(id_node);
	assert(pnode != NULL);

	// intilization
	dForce = 0.0;
	num = 0;

	// right hand
	x = (int)floor(pnode->GetX());
	y = (int)floor(pnode->GetY());

	for (j=0; j<=0; j++)
	{
		x1 = j+x;
		for (i=0; i<=0; i++)
		{
			y1 = i+y;
			if (x1>=0 && x1<=player_region->height-1 && y1>=0 && y1<=player_region->width)
				num++;
		}
	}
	if (num>0)
		dForce = GetLogPostProb(x,y);
	else
		dForce = 100.0;
	return dForce;
}


McPoint2D GetOffsetFromNormalAndForces(const double normal_angle,	const double sfs, const double curvature,
									 double* pdStatForce, double* pdSmoothForce,
									 McPoint2D* pofStatForce, McPoint2D* pofSmoothForce)
{
	double		dDf=FORCE_NORMAL,fs;
	McPoint2D	pt2dOffset;

	// initialization
	if (pofStatForce != NULL)
		pofStatForce->Set(0.0,0.0);
	if (pofSmoothForce != NULL)
		pofSmoothForce->Set(0.0,0.0);

	fs = (-sfs+0.5*MU*curvature)/dDf;

	pt2dOffset.Set(fs*cos(normal_angle), fs*sin(normal_angle));

	if (pdStatForce != NULL)
		*pdStatForce = -sfs/dDf;
	if (pdSmoothForce != NULL)
		*pdSmoothForce = -0.5*MU*curvature/dDf;

	if (pofStatForce != NULL)
	{
		double dStatForce;

		dStatForce = -sfs/dDf;
		pofStatForce->Set(dStatForce*cos(normal_angle),dStatForce*sin(normal_angle));
	}

	if (pofSmoothForce != NULL)
	{
		double dSmoothForce;
		dSmoothForce = +0.5*MU*curvature/dDf;
		pofSmoothForce->Set(dSmoothForce*cos(normal_angle),dSmoothForce*sin(normal_angle));
	}

	return pt2dOffset;
}


McPoint2D MCMC_AtomicRegion::GetDiffusionMove(const int id_node1, const int id_node, const int id_node2,
											  McPoint2D &ofStatForce, McPoint2D &ofSmoothForce)
{
	double		dnorm,curv,sfs,dmin_dis;
	McPoint2D	pt2dOffset;

	// initialization
	ofStatForce.Set(0.0,0.0);
	ofSmoothForce.Set(0.0,0.0);

	if (id_node1>=0 && id_node>=0 && id_node2>=0)
	{
		// get norm direction and curvature
		dnorm = GetNormal(id_node1,id_node,id_node2,curv,&dmin_dis);
		// get statistics force
		sfs = GetStatForce(id_node);

		pt2dOffset = GetOffsetFromNormalAndForces(dnorm, sfs, curv, NULL, NULL, &ofStatForce, &ofSmoothForce);
	}
	else
		// didn't find any vertex
		pt2dOffset.Set(0.0,0.0);

	return pt2dOffset;
}

McPoint2D MCMC_AtomicRegion::GetDiffusionMove(const int index, McPoint2D &ofStatForce, McPoint2D &ofSmoothForce)
{
	int id_node1,id_node,id_node2;

	// make sure id_node_cl1, id_node and id_node2_cl appear in clockwise order
	id_node1 = GetSortedLastNode(index);
	id_node = GetSortedNode(index);
	id_node2 = GetSortedNextNode(index);

	return GetDiffusionMove(id_node1,id_node,id_node2,ofStatForce,ofSmoothForce);
}


int	MCMC_AtomicRegion::GetCorrespondingVertex(const int id_node)
{
	int id_cseg,id_ver;

	id_ver = -1;
	id_cseg = GetCurrentCSeg();
	if (id_cseg>=0)
	{
		MCMC_Contour_Segment *pcseg;
		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);
		if (id_node==pcseg->node[0])
			id_ver = pcseg->end[0];
		else if (id_node==pcseg->node[1])
			id_ver = pcseg->end[1];
	}
	return id_ver;
}

// get the bound for this region
irmm_Rect MCMC_AtomicRegion::GetBound(irmm_Rect *pbound)
{
	irmm_Rect	bound;
	int			id_node;
	MCMC_Node	*pnode;

	if (pbound==NULL)
		bound.Set(9999,9999,0,0);
	else
		bound = *pbound;

	id_node = GetFirstNode();
	while (id_node>=0)
	{
		pnode = LayerRegion()->GetNode(id_node);
		bound.SetBound((int)pnode->GetY()-2,(int)pnode->GetX()-2);
		bound.SetBound((int)pnode->GetY()+2,(int)pnode->GetX()+2);

		id_node = GetNextNode();
	}
	if (bound.left<0)					bound.left = 0;
	if (bound.top<0)					bound.top = 0;
	if (bound.right>=LayerRegion()->width)	bound.right = LayerRegion()->width-1;
	if (bound.bottom>=LayerRegion()->height)	bound.bottom = LayerRegion()->height-1;

	return bound;
}


void MCMC_AtomicRegion::UpdateTopologyRecord(const bool bupdate_label_map)
{
	MCMC_Contour_Segment *pcseg;

	pcseg = player_region->GetCSeg(GetFirstCSeg());
	while (pcseg != NULL)
	{
		pcseg->UpdateTopologyRecord();
		pcseg = player_region->GetCSeg(GetNextCSeg());
	}

	// label map
	if (bupdate_label_map)
		Label(player_region->mx_label_id_aregions, Id());
}

void MCMC_AtomicRegion::GetAllAtomicRegions(LinkList<int> &list_aregions)
{
	list_aregions.Empty();
	list_aregions.Insert(Id());
}


void MCMC_AtomicRegion::UpdateLowLevelStat(MCMC_FilterResp &resp)
{
	assert(plow_level_stat != NULL);
	plow_level_stat->Update(resp, pixel_set);
}


///////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

// construction & destruction
MCMC_GRegion::MCMC_GRegion()
{
}

MCMC_GRegion::~MCMC_GRegion()
{
	Empty();
}

void MCMC_GRegion::EmptyTopology(void)
{
	list_indexid.Empty();
	list_cseg_start.Empty();

	btopology_consistent = true;
	list_all_cseg.Empty();
}

void MCMC_GRegion::Empty(void)
{
	MCMC_RegionElement::Empty();
	
	EmptyTopology();
}

MCMC_GRegion& MCMC_GRegion::operator=(MCMC_GRegion& a)
{
	if (this != &a)
	{
		*((MCMC_RegionElement*) this) = *((MCMC_RegionElement*) &a);

		list_indexid	= a.list_indexid;
		list_cseg_start	= a.list_cseg_start;		        //        ___ 
		                                                    //       |___|
		// for g-region with dangling atomic region               ---
		btopology_consistent = a.btopology_consistent;  //   |   |
		list_all_cseg = a.list_all_cseg;                    //    ---
	}

	return *this;
}


ofstream &	operator>>(MCMC_GRegion &a, ofstream &strm)
{
	*((MCMC_RegionElement*)&a) >> strm;
	*((MCMC_MemoryNode*)&a) >> strm;

	a.list_indexid >> strm;
	a.list_cseg_start >> strm;

	return strm;
}

ifstream &	operator<<(MCMC_GRegion &a, ifstream &strm)
{
	a.Empty();

	*((MCMC_RegionElement*)&a) << strm;
	*((MCMC_RegionElement*)&a) << strm;

	a.list_indexid << strm;
	a.list_cseg_start << strm;

	return strm;
}

int	MCMC_GRegion::GetFirstVer(void)
{
	if (btopology_consistent)
		return MCMC_RegionElement::GetFirstVer();
	else
	{
		if (!Valid())
			return -1;
		
		int	current_seg=-1;

		current_seg = -1;
		list_all_cseg.GetFirst(current_seg);

		MCMC_Contour_Segment* pcseg = LayerRegion()->GetCSeg(current_seg);
		assert(pcseg != NULL);
		if (pcseg==NULL)
			return -1;
		int iVer1,iVer2;
		pcseg->GetVertices(iVer1,iVer2);
		if (IsTheSameElement(pcseg->GetRightAtomicRegion(), IndexId()))
			current_ver = iVer1;
		else
			current_ver = iVer2;

		last_seg = -1;
		return	current_ver;
	}
}

int	MCMC_GRegion::GetNextVer(void)
{
	if (btopology_consistent)
		return MCMC_RegionElement::GetNextVer();
	else
	{
		if (!Valid())
			return -1;

		int	current_seg;

		current_seg = -1;
		list_all_cseg.GetNext(current_seg);

		if (current_seg<0)
			return -1;
		last_seg = current_seg;
		// get corresponding contour segment
		MCMC_Contour_Segment* pcseg = LayerRegion()->GetCSeg(current_seg);
		assert(pcseg != NULL);
		// go to the next vertex
		current_ver = pcseg->GetTheOtherVertex(current_ver);

		if (current_seg==start_cseg || current_seg<0)	// go back to starting point
			return -1;
		else
			return current_ver;
	}
}

int	MCMC_GRegion::GetCurrentCSeg(void)
{
	if (IsTopologyConsistent())
		return MCMC_RegionElement::GetCurrentCSeg();
	else
	{
		int current_seg=-1;

		if (current_ver >= 0)
			list_all_cseg.GetCurrent(current_seg);
		
		return current_seg;
	}
}

int	MCMC_GRegion::GetNextNode(void)
{
	if (IsTopologyConsistent())
		return MCMC_RegionElement::GetNextNode();


	// get current segment
	int id_cseg = GetCurrentCSeg();
	if (id_cseg < 0)
		return -1;
	else
	{
		int					 id_node;
		MCMC_Contour_Segment *pcseg;

		pcseg = LayerRegion()->GetCSeg(id_cseg);
		assert(pcseg != NULL);

		if (direction)
			id_node = pcseg->GetNextNode();
		else
			id_node = pcseg->GetPrevNode();
		if (id_node >= 0)
		{
			id_node_temp = id_node;
			if (id_node != id_node_first)
				return id_node;		// get the next node
			else
				return -1;
		}
		else
		{
			// start from the next segment
			id_cseg = GetNextCSeg();
			if (id_cseg < 0)
				return -1;
			else
			{
				pcseg = LayerRegion()->GetCSeg(id_cseg);
				assert(pcseg != NULL);

				// skip the starting node
				id_node = pcseg->GetFirstNode();
				if (id_node == id_node_temp)
					direction = 1;
				else
				{
					id_node = pcseg->GetLastNode();
					direction = 0;
				}
				if (id_node == id_node_temp)
				{
					if (direction)
						id_node = pcseg->GetNextNode();
					else
						id_node = pcseg->GetPrevNode();
				}
				else
					id_node = pcseg->GetFirstNode();

				id_node_temp = id_node;
				if (id_node != id_node_first)
					return id_node;
				else
				{
					if (list_all_cseg.IsAtTail())
						return -1;
					else
						return id_node;
				}
			}
		}
	}
}


void MCMC_GRegion::GetCSetFromVertex(MCMC_Vertex *pver, mcmcVector<int> &vc_csegs, const int id_region, const int side)
{
	if (pver != NULL)
		pver->GetCSegOfGRegion(vc_csegs, id_region, side);
}

bool MCMC_GRegion::IsTheSameElement(const int id_aregion1, const int indexid_region2)
{
	return LayerRegion()->IsInGRegion(id_aregion1, indexid_region2);
}

// depth first search
int	MCMC_GRegion::GetFirstAtomicRegion(void)
{
	bool				bgo;
	int					indexid_region,id_aregion;
	MCMC_RegionElement	*pregion_element;

	bgo = list_indexid.GetFirst(indexid_region);
	while (bgo)
	{
		pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_region);
		assert(pregion_element != NULL);

		id_aregion = pregion_element->GetFirstAtomicRegion();
		if (id_aregion >= 0)
			return id_aregion;

		bgo = list_indexid.GetNext(indexid_region);
	}
	return -1;
}

int	MCMC_GRegion::GetCurrentAtomicRegion(void)
{
	bool				bgo;
	int					indexid_region,id_aregion;
	MCMC_RegionElement	*pregion_element;

	bgo = list_indexid.GetCurrent(indexid_region);
	if (!bgo)
		return -1;

	pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_region);
	assert(pregion_element != NULL);

	id_aregion = pregion_element->GetCurrentAtomicRegion();
	if (id_aregion >= 0)
		return id_aregion;
	else
		return -1;
}


// depth first search
int MCMC_GRegion::GetNextAtomicRegion(void)
{
	bool				bgo;
	int					indexid_region,id_aregion;
	MCMC_RegionElement	*pregion_element;

	bgo = list_indexid.GetCurrent(indexid_region);
	if (!bgo)
		return -1;

	pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_region);
	assert(pregion_element != NULL);

	id_aregion = pregion_element->GetNextAtomicRegion();
	if (id_aregion >= 0)
		return id_aregion;

	bgo = list_indexid.GetNext(indexid_region);
	while (bgo)
	{
		pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_region);
		assert(pregion_element != NULL);

		id_aregion = pregion_element->GetFirstAtomicRegion();
		if (id_aregion >= 0)
			return id_aregion;

		bgo = list_indexid.GetNext(indexid_region);
	}

	return -1;
}

int	MCMC_GRegion::NumAtomicRegions(void)
{
	int					indexid_region,num;
	bool				bgo;
	MCMC_RegionElement	*pregion_element;

	num = 0;
	bgo = list_indexid.GetFirst(indexid_region);
	while (bgo)
	{
		pregion_element = LayerRegion()->GetRegionFromIndexId(indexid_region);
		assert(pregion_element != NULL);

		num += pregion_element->NumAtomicRegions();

		bgo = list_indexid.GetNext(indexid_region);
	}
	return num;
}


void MCMC_GRegion::FindTheStartCSeg(void)
{
	// all possible starting contour segments
	int num;
	LinkList<int> *p_lists = FindClosedContours(num);
	if (p_lists != NULL)
		delete []p_lists;
}

bool MCMC_GRegion::UseFirstStartCSeg(void)
{
	if (list_cseg_start.Num() <= 1)
		return true;
	else
	{
		int temp;

		list_cseg_start.GetFirst(temp);
		start_cseg = temp;
		return true;
	}
}

bool MCMC_GRegion::UseNextStartCSeg(void)
{
	if (list_cseg_start.Num() <= 1)
		return false;
	else
	{
		int	temp;

		if (list_cseg_start.GetNext(temp))
		{
			start_cseg = temp;
			return true;
		}
		else
		{
			list_cseg_start.GetFirst(temp);
			start_cseg = temp;
			return false;
		}
	}
}

int	MCMC_GRegion::GetFirstRegion(void)
{
	int	temp;
	if (list_indexid.GetFirst(temp))
		return temp;
	else
		return -1;
}

int	MCMC_GRegion::GetNextRegion(void)
{
	int	temp;
	if (list_indexid.GetNext(temp))
		return temp;
	else
		return -1;
}

bool MCMC_GRegion::GetFirstPixel(int &j, int &i)
{
	MCMC_AtomicRegion *paregion;

	// retrieve the atomic region
	paregion = LayerRegion()->GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		if (paregion->GetFirstPixel(j, i))
			return true;

		paregion = LayerRegion()->GetAtomicRegion(GetNextAtomicRegion());
	}

	return false;
}

bool MCMC_GRegion::GetNextPixel(int &j, int &i)
{
	MCMC_AtomicRegion *paregion;

	// retrieve the atomic region
	paregion = LayerRegion()->GetAtomicRegion(GetCurrentAtomicRegion());
	if (paregion->GetNextPixel(j, i))
		return true;

	paregion = LayerRegion()->GetAtomicRegion(GetNextAtomicRegion());
	while (paregion != NULL)
	{
		if (paregion->GetFirstPixel(j, i))
			return true;

		paregion = LayerRegion()->GetAtomicRegion(GetNextAtomicRegion());
	}
	return false;
}

void MCMC_GRegion::GetAllAtomicRegions(LinkList<int> &list_aregions)
{
	list_aregions.Empty();

	int id_aregion;

	id_aregion = GetFirstAtomicRegion();
	while (id_aregion >= 0)
	{
		list_aregions.Insert(id_aregion);
		id_aregion = GetNextAtomicRegion();
	}
}

// create the pxiel set for this region
void MCMC_GRegion::CreateRegionDomain(void)
{
	MCMC_AtomicRegion *paregion;
	bool	bgo;
	int		i,j;
	int	num_pixels;

	// initialization
	pt_center.Set(0.0, 0.0);
	num_pixels = 0;

	// retrieve each region
	paregion = LayerRegion()->GetAtomicRegion(GetFirstAtomicRegion());
	while (paregion != NULL)
	{
		paregion->CreateRegionDomain();
		
		bgo = paregion->GetFirstPixel(j, i);
		while (bgo)
		{
			pt_center.x += j;
			pt_center.y += i;
			num_pixels++;

			bgo = paregion->GetNextPixel(j, i);
		}
		paregion = LayerRegion()->GetAtomicRegion(GetNextAtomicRegion());
	}

	if (num_pixels > 0)
		pt_center = pt_center/((double)num_pixels);
	
}

void MCMC_GRegion::AddRegion(const int indexid)
{
	MCMC_RegionElement *pregion;
	pregion = LayerRegion()->GetRegionFromIndexId(indexid);
	if (pregion == NULL)
		return;

	int temp=indexid;
	
	list_indexid.Insert(temp);
	pregion->SetIndexIdParent(IndexId());
}

void MCMC_GRegion::DeleteRegion(const int indexid)
{
	int temp=indexid;
	
	list_indexid.Delete(temp);

	MCMC_RegionElement *pregion;
	pregion = LayerRegion()->GetRegionFromIndexId(indexid);
	if (pregion != NULL)
		pregion->SetIndexIdParent(-1);
}

// get the length of its total boundary
double MCMC_GRegion::GetLength(void)
{
	MCMC_Contour_Segment *pcseg;
	int		id_cseg,id_cseg_buf,id_temp;
	double	dlength=0.0;
	bool	bgo;
	int		max_cseg,k;

	
	max_cseg = LayerRegion()->NumCSegments();

	id_cseg_buf = start_cseg;
	
	k = 0;
	bgo = list_cseg_start.GetFirst(id_temp);
	while (bgo)
	{
		start_cseg = id_temp;
		id_cseg = GetFirstCSeg();
		while (id_cseg>=0)
		{
			pcseg = LayerRegion()->GetCSeg(id_cseg);
			assert(pcseg != NULL);

			dlength += pcseg->GetLength();
			id_cseg = GetNextCSeg();
			k++;

			if (k > max_cseg)
			{
				cout << "GetLength" << endl;
				XPT_OutOfBound *pxpt;
				pxpt = new XPT_OutOfBound("Wrong topology.");				
				throw pxpt;
			}
		}
		// change the starting contour segment
		bgo = list_cseg_start.GetNext(id_temp);
	}
	
	
	start_cseg = id_cseg_buf;
	return dlength;
}


LinkList<int>* MCMC_GRegion::FindClosedContours(int &num)
{
	LinkList<int> list_aregions;
	
	// initialization
	btopology_consistent = true;
	list_all_cseg.Empty();
	
	
	GetAllAtomicRegions(list_aregions);

	// get the contour segments on the boundary
	LinkList<int>			list_cseg_boundary;
	MCMC_Contour_Segment	*pcseg;
	MCMC_AtomicRegion		*paregion,*paregion_left,*paregion_right;
	bool					bgo;
	int						id_aregion;

	bgo = list_aregions.GetFirst(id_aregion);
	while (bgo)
	{
		paregion = LayerRegion()->GetAtomicRegion(id_aregion);
		
		// retrieve each contour segment
		pcseg = LayerRegion()->GetCSeg(paregion->GetFirstCSeg());
		while (pcseg != NULL)
		{
			paregion_left = LayerRegion()->GetAtomicRegion(pcseg->GetLeftAtomicRegion());
			paregion_right = LayerRegion()->GetAtomicRegion(pcseg->GetRightAtomicRegion());
			
			if (paregion_left==NULL || paregion_right==NULL)
				list_cseg_boundary.Insert(pcseg->Id());
			else
			{
				if (!IsTheSameElement(paregion_left->Id(), IndexId()) || 
					!IsTheSameElement(paregion_right->Id(), IndexId()))
					list_cseg_boundary.Insert(pcseg->Id());
			}
			pcseg = LayerRegion()->GetCSeg(paregion->GetNextCSeg());
		}
		bgo = list_aregions.GetNext(id_aregion);
	}
	

	list_cseg_start.Empty();
	
	int				id_cseg,num_contours,icseg_longest;
	LinkList<int>	list_temp;
	double	dlength_temp,dlength_longest;
	bool bconsistent=true;

	// initialization
	num_contours = 0;
	icseg_longest = 0;
	
	dlength_temp = 0.0;
	dlength_longest	 = -1.0;

	// obtain each closed contour
	while (list_cseg_boundary.Num() > 0)
	{
		list_cseg_boundary.GetFirst(id_cseg);
		ChangeStartCSeg(id_cseg);

		dlength_temp = 0.0;

		id_cseg = GetFirstCSeg();
		while (id_cseg >= 0)
		{
			pcseg = player_region->GetCSeg(id_cseg);
			assert(pcseg != NULL);
			dlength_temp += pcseg->GetLength();

			if (!list_cseg_boundary.Delete(id_cseg))
			{
				id_cseg = -1;
				bconsistent =false;
				break;
			}
			else
				list_temp.Insert(id_cseg);

			id_cseg = GetNextCSeg();
		}
		// Is this contour the longest?
		if (dlength_temp > dlength_longest)
		{
			icseg_longest = num_contours;
			dlength_longest = dlength_temp;
		}

		list_temp.Insert(id_cseg);		
		num_contours++;
	}
	
	btopology_consistent = bconsistent;
	num = 0;
	if (num_contours == 0)
		return NULL;
	else
	{
		LinkList<int> *p_lists;
		bool bfirst=true;

		p_lists = new LinkList<int>[num_contours];
		bgo = list_temp.GetFirst(id_cseg);
		while (bgo)
		{
			if (id_cseg < 0)
			{
				num++;
				bfirst = true;
			}
			else
			{
				p_lists[num].Insert(id_cseg);
				if (bfirst)
				{
					list_cseg_start.Insert(id_cseg);
					if (num == icseg_longest)
						ChangeStartCSeg(id_cseg);	// setup this one as the starting cseg

					bfirst = false;
				}
			}
			bgo = list_temp.GetNext(id_cseg);
		}
		if (icseg_longest != 0)
		{
			assert(list_cseg_start.Delete(start_cseg));
			list_cseg_start.InsertToHead(start_cseg);

			list_temp = p_lists[0];
			p_lists[0] = p_lists[icseg_longest];
			p_lists[icseg_longest] = list_temp;
		}
		
		// copy all the contour segments into the list
		if (!btopology_consistent)
		{
			for (int i=0; i<num; i++)
				list_all_cseg = list_all_cseg+p_lists[i];

			int temp=-1;
			list_all_cseg.GetFirst(temp);
			start_cseg = temp;

			list_cseg_start.Empty();
			list_cseg_start.Insert(start_cseg);
		}
		return p_lists;
	}
}

void MCMC_GRegion::GetMinDistTo(mcmcVector<double> &vc_dist, SimpleMatrix<McPoint2D> &vc_points)
{
	// to compute the minimal distance
	LinkList<double> list;
	double dist;
	int m=0;

	for (int i=0; i<vc_points.size()/2; i++)
	{
		dist = MinDistToPoints(vc_points, i);
		list.Insert(dist);
		if (dist > 6.0)
			m++;
		if (m > 5)
			break;
	}

	list.ToVector(vc_dist);
}

double MCMC_GRegion::ComputeEnergyShape(void)
{
	energy_shapeprior = GetLength()*PRIOR_SHAPE + PRIOR_REGION;

	// penalty for disjoint components
	if (list_cseg_start.Num() > 1)
		energy_shapeprior += PRIOR_REGION*3.5*(list_cseg_start.Num()-1);
	if (!IsTopologyConsistent())
		energy_shapeprior += PRIOR_REGION*2.0;

	bool				bgo,bgo2;
	LinkList<McPoint2D>	list_points;
	MCMC_Node			*pnode;
	McPoint2D			pt_node;
	SimpleMatrix<McPoint2D>	vc_points;
	mcmcVector<double>		vc_dist;
	int					i,num_not_to_consider,k,m;
	double				dmean,dvar,dmean_robust;

	
	bgo = UseFirstStartCSeg();
	while (bgo)
	{
		list_points.Empty();
		
		// obtain all the points
		pnode = LayerRegion()->GetNode(GetFirstNode());
		while (pnode != NULL)
		{
			pt_node.Set(pnode->GetX(), pnode->GetY());
			list_points.Insert(pt_node);

			pnode = LayerRegion()->GetNode(GetNextNode());
		}
	
		if (list_points.Num() > 15 && list_points.Num()<2*(max(LayerRegion()->height, LayerRegion()->width)+5))
		{
			LinkList<McPoint2D>	list_sampled_points;
			if (list_points.Num()>80)
			{
				k = (int)floor(((double)list_points.Num())/60);
				m = 0;
				bgo2 = list_points.GetFirst(pt_node);
				while (bgo2)
				{
					if (m%k == 0)
						list_sampled_points.Insert(pt_node);
					m++;
					bgo2 = list_points.GetNext(pt_node);
				}

				list_sampled_points.ToSimpleMatrix(vc_points);
			}
			else
				list_points.ToSimpleMatrix(vc_points);

			dArea = area(vc_points);

			if (Area()/(vc_points.size()/2) < 1/6.0*vc_points.size()) // elongated
			{
				GetMinDistTo(vc_dist, vc_points);
				dmean = vc_dist.Mean();

				if (dmean < MAX_DIST_TO_CONSIDER/3.0*2.0)
				{
					// points may be in the extreme
					num_not_to_consider = min(5, IntClose(dmean));

					k = 0;
					m = 0;
					dmean_robust = 0.0;
					for (i=0; i<vc_dist.size(); i++)
					{
						if (vc_dist(i) < dmean+5.0 || m>=num_not_to_consider)
						{
							dmean_robust += vc_dist(i);
							k++;
						}
						else
							m++;
					}
					if (k > 0)
					{
						dmean_robust /= k;
						dvar = 0.0;
						m = 0;
						for (i=0; i<vc_dist.size(); i++)
						{
							if (vc_dist(i) < dmean+5.0 || m>=num_not_to_consider)
							{
								dvar += pow(vc_dist(i)-dmean_robust, 2.0);
								k++;
							}
							else
								m++;
						}
						dvar = sqrt(dvar)/k;
						
						if (dmean_robust < 5.5 && dvar < 0.17)
							energy_shapeprior -= PRIOR_REGION*3.0;
					}
				}
			}
		}
		bgo = UseNextStartCSeg();
	}

	return energy_shapeprior;
}

double MCMC_GRegion::ComputeEnergyShapeForList(LinkList<int> &list_aregions)
{
	if (list_aregions.Num() <= 0)
		return 0.0;
	
	LinkList<int>	list_record_idparent;
	int				id_aregion,indexid_parent;
	bool			bgo;
	MCMC_AtomicRegion *paregion;
	
	// for the topology part
	EmptyTopology();

	// save the parent id of these atomic regions
	bgo = list_aregions.GetFirst(id_aregion);
	while (bgo)
	{
		paregion = LayerRegion()->GetAtomicRegion(id_aregion);
		assert(paregion != NULL);

		indexid_parent = paregion->IndexIdParent();
		list_record_idparent.Insert(indexid_parent);

		AddRegion(paregion->IndexId());
			
		bgo = list_aregions.GetNext(id_aregion);
	}


	LinkList<int>	*p_lists;
	int		num;

	p_lists = FindClosedContours(num);

	// restore
	assert(list_aregions.Num() == list_record_idparent.Num());

	bgo = list_aregions.GetFirst(id_aregion);
	list_record_idparent.GetFirst(indexid_parent);
	while (bgo)
	{
		paregion = LayerRegion()->GetAtomicRegion(id_aregion);
		assert(paregion != NULL);
		paregion->SetIndexIdParent(indexid_parent);

		bgo = list_aregions.GetNext(id_aregion);
		list_record_idparent.GetNext(indexid_parent);
	}
	
	if (p_lists != NULL)
		delete []p_lists;
	return ComputeEnergyShape();
}


void MCMC_GRegion::SetParentIdForChildRegions(void)
{
	MCMC_RegionElement *pregion;
	
	pregion = LayerRegion()->GetRegionFromIndexId(GetFirstRegion());
	while (pregion != NULL)
	{
		pregion->SetIndexIdParent(IndexId());
		pregion = LayerRegion()->GetRegionFromIndexId(GetNextRegion());
	}
}

void MCMC_GRegion::UpdateLowLevelStat(MCMC_FilterResp &resp)
{
}
